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1. Introduction 


Vision-based ranging has emerged as the central problem in automating the helicopter nap-of- 
the-earth flight regime. Current status of the research in vision-based ranging at NASA Ames is 
summarized in Reference 1. The present work is a component of this research effort. 

The recovery of three dimensional geometry from two dimensional scenes is based on the fact that 
if one has at least two views of a scene obtained from different vantage points, the difference between 
the location of corresponding objects in the images is a measure of the range. The process of finding 
the corresponding objects within the field of view is the correspondence problem in machine vision 
literature. Correspondence problem can be solved if all the images unambiguously contained every 
object in the scene. Since real imaging devices have limited field of view and resolution, concept of 
correspondence is in reality a correspondence hypothesis. The relative object displacement obtained by 
satisfying the correspondence hypothesis in these images are called disparities. The farther the objects 
are with respect to the imaging device, the less the observed disparity would be. For instance, objects 
that are very far away would appear at nearly the same location in the images, while closer objects 
will exhibit a large disparity. If the disparity of the corresponding objects in an image sequence could 
be measured, the perspective projection equations could be used to compute the scene depth. This 
fact forms the basis for a large majority of vision-based ranging methods discussed in the literature. 
Sometimes, establishing correspondence between several views of the scene may improve the range 
estimation accuracy. 

Two distinct families of techniques can be identified for satisfying the correspondence hypothesis. 
Methods which establish correspondence between relativly few chosen features of interest in the images 
are termed feature based methods. On the other hand, techniques which attempt to establish the 
correspondence between every point in the images are called field-based methods. The focus of present 
research effort is on field-based ranging technique. The present research has its basis in the robot vision 
literature [2]. 

In the present work, a family of field-based ranging algorithms are cast in an analytical format. 
An advantage of the analytical approach is that it reveals the means for systematic performance 
improvement. If desired, additional heuristics may be included in the formulation to enhance the 
algorithm performance. 

These ranging algorithms have their basis in a mathematical approximation of the correspondence 
hypothesis together with the perspective projection geometry. The central ideas in the present work 
are: 

• Mathematical approximation of the correspondence hypothesis 

• The use of incremental perspective projection in the correspondence hypothesis approximation 

to yield vision-based ranging equations 

• Methods for estimating ranging equation parameters 

Ensuing chapters in this report are based on a few papers written over the past three years. The 
motivation and summary of the research is given at the beginning of each chapter, followed by the 
paper. The contributions of the present research are given in a condensed form in the Conclusions 
section. 

A part of the research effort was expended in synthesizing a guidance law for performing obstacle 
avoidance using the image-based sensor data. The range data generated by image-based sensors are 
generally discrete and are available only at about 10 % of the points on the image. An optimal 
guidance law that explicitly uses the discrete range data together with a nonlinear point-mass vehicle 
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model was formulated using the calculus of variations. Details of this guidance law development and 
implementation are given in Chapter 6. 


2. Development of a Field-Based Ranging Equation 

Initial phases of the present research focussed on ranging using motion image sequences generated 
by a single camera fixed to a moving vehicle. A ranging scheme using temporal image sequences is 
developed based on the Optical Flow Constraint Equation of Horn and Schunck [3]. The optical flow 
constraint equation relates the temporal partial derivatives with the spatial partial derivatives of the 
image irradiance E as: 


dE 

dt 


dE dE 
+ ^—u + —v = 0 

OX„ OZr, 


( 2 . 1 ) 


where (x p , z p ) are the major and minor axes of the image plane while u and v are called the optical 
flow components. Generally, the ranging problem is treated in a two-step fashion. First, the optical 
flow components u and v are computed. Since just one expression is available for the calculation of 
two quantities, an additional condition has to be imposed to obtain a unique solution. Horn and 
Schunk [3] propose the use of a smoothness constraint for the optical flow components to alleviate this 
difficulty. The optical flow components are then used for computing the range. 

Instead of this two step process, incremental perspective projection can be combined directly with 
the optical flow constraint to yield a ranging equation. Since such a concatenation is direct for the 
case of purely translational motion of the camera, this case was investigated first. 

The perspective projection relationships for a camera located at a position (xq, J/o) are given by 


_ f(x - x 0 ) _ f{z - z 0 ) 

p (y-yo) ’ 2p " (y-yo) 

These equations can be differentiated to yield the optical flow components as 


( 2 . 2 ) 


u = 


[x p yo - fi 0 ] 


V — 


{Zpijo - fz o] 


(2.3) 


( y ~ yo) “ (y - yo) 

Substituting these in the optical flow constraint equation yields an Optical Navigation Equation 


as : 


dE 

dt 


^~~( x pyo ~ A o) + irr( z pyo ~ Ao) 


dz-r, 


(y - yo) 


(2.4) 


This equation can be used to compute the range (y — y Q ) from the image irradiance partial deriva- 
tives and camera position components. A simulated scene was used to test the algorithm performance. 
A paper based on this research was presented at the 1989 AIAA Guidance, Navigation, and Control 
Conference discussing the details of this investigation. This paper is included in this chapter. 
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Passive Navigation Using Image Irradiance Tracking 


P. K. A. Menon'and B. Sridhar* 
NASA Ames Research Center 
Moffett Field, CA 94035 


Abstract 

Rotorcraft operating at low altitudes require nav- 
igational schemes for locating the terrain and ob- 
stacles. Due to the covert nature of missions to be 
accomplished, a passive navigation scheme is de- 
sirable. This paper describes the development of 
a passive navigation scheme combining image se- 
quences from a vehicle mounted camera with ve- 
hicle motion variables. Geometric properties of 
perspective projection together with an image ir- 
radiance tracking scheme at each pixel are used 
to determine the range to various objects within 
the field-of-view. Derivation of the numerical al- 
gorithm and simulation results are given. Other 
applications of the proposed approach include nav- 
igation for autonomous planetary rovers and teler- 
obots. 

1 Introduction 

The problem of determining the scene geome- 
try from on-board sensors is emerging as a cen- 
tral issue in nap-of-the-earth helicopter flight 
guidance 1 ’ 2 . The objective of passive naviga- 
tion is to use image-based sensors such as FLIR 
and low-light level TV to determine the range 
to various objects within the field-of-view. Sev- 
eral approaches to the passive navigation problem 
have been suggested in recent literature 1,3-6 . Ap- 
proaches for passive navigation using stereo image 
pairs have also been reported 7 . However, with the 
exception of Reference 1, none of these techniques 

‘Member AIAA, Georgia Institute of technology, At- 
lanta. Mailing Address : FSN Branch, MS 210-9 

Research Scientist, Member AIAA, FSN Branch, MS 
210-9 


appear to have been developed to a point where 
they can be used for vehicle navigation. Moreover, 
there exists an ambiguity regarding the interpre- 
tation of the visual navigation task. With the ex- 
ception of Reference 1, previous researchers have 
tended to interpret image based passive naviga- 
tion as the process of simultaneously determining 
the vehicle motion find object position parameters 
from an image sequence. This can be shown to 
produce an under- determined system and the so- 
lution can be obtained only by imposing an addi- 
tional super-criterion. The present approach, how- 
ever, assumes that the vehicle motion parameters 
are known from the on-board inertial navigation 
system. As a result, the navigation task is to use 
the available image sequence together with the in- 
ertial navigational information to determine the 
relative location of various objects within the field 
of view. 

The data available in a monochromatic image 
is a distribution of image irradiance specified as 
a function of imaging device aperture coordinates. 
Various points in an image obtained from a mov- 
ing vehicle will exhibit an irradiance distribution 
change as a function of time. This change depends 
on the relative location of various objects, vehicle 
motion parameters and the surface reflectance. If 
the surface reflectance is assumed to remain con- 
stant, then the observed image irradiance distri- 
bution change is entirely due to the relative loca- 
tion of various objects within the field-of-view and 
vehicle motion parameters. In this case, if the ve- 
hicle motion parameters and the imaging device 
constants are known, it is possible to determine 
the realtive location of various points within the 
field of view. At this early stage in research, the 
investigation will be restricted to objects that are 
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stationary with respect to the earth fixed frame. 
Additionally, only pure translational motion of the 
vehicle will considered. 

The approach presented in this paper has its 
roots in the optical Sow constraint equation de- 
rived by Horn find Schunck 3,4 ’ 8 . With the as- 
sumption that the objects within the field-of-view 
are fixed, this equation is combined with perspec- 
tive projection geometry and the inertial naviga- 
tion system data to derive a direct navigational 
equation. This navigational equation relates the 
spatial and temporal partial derivatives of the im- 
age irradiance with the distance to various objects 
within the field-of-view and the vehicle motion pa- 
rameters. An expression for determining the spa- 
tial and temporal sampling intervals required for 
satisfactory performance of the proposed scheme 
is next derived. Details of this analysis are given 
in Section 2. In Section 3, an optimal approach 
for estimating various partial derivatives is syn- 
thesized using the Calculus of Variations 9,10 . This 
approach requires the solution of a linear, second- 
order partial differential equation. Using the exist- 
ing theory of such partial differential equations 11 , 
the solution is obtained by assuming a product 
decomposition. Spatial discretization of this so- 
lution then produces feedback scheme for track- 
ing the image irradiance as well as a systematic 
approach for estimating the spatial and temporal 
partial derivatives. 

The performance of the proposed navigation 
scheme is demonstrated using a simulated image 
sequence in Section 4. Conclusions and future re- 
search directions are indicated in Section 5. 


X p — Z p plane, with E being the irradiance speci- 
fied on a gray scale at a point (x p , z p ) on the image. 
To avoid the complications that can arise while 
considering an inverted image as is encountered in 
real cameras, the front image plane 12 will be used 
for image description. It is assumed in the ensu- 
ing analysis that the image irradiance distribution 
is continuous except at finite number of regions. 
Thus, for the purposes of the present analysis, ef- 
fects such as spatial image discretization and gray 
scale quantization are ignored. The image plane is 
assumed to be located at a distance from 

the origin of the vehicle body axis system, and 
oriented with respect to the body axis by three 
rotations €i,€ 2 , The vehicle body axis system 
may be related to an earth fixed frame through 
three Euler angles Vs 9A and three translations 
yo> ^o. The vehicle inertial navigation system is 
assumed to produce estimates of a^yo*^ t/>, 0,<£ 
and their time derivatives. Let the location of 
the various objects within the field-of-view with 
respect to the earth fixed frame be x, y, z. Various 
coordinate systems outlined in the foregoing are 
illustrated in Figure 1, 

Any point sc p ,z p corresponding to an object 
in the image plane is related to its location with 
respect to the earth-fixed frame through various 
translations and rotations, together with the per- 
spective projection. Considering pure transla- 
tions, if one assumes that the camera is located 
at the origin of the vehicle body coordinate sys- 
tem and aligned with it, the relation between an 
image point and its position with respect to the 
earth fixed frame is given by 


2 The Navigation Scheme 

Assuming that the image is rectangular with the 
major axis parallel to the horizontal, one may de- 
fine the major axis as the X p axis and the minor 
axis as the Z p axis. The Y p axis is normal to the 
image plane. The origin of this image plane coor- 
dinate system may be at the center of the aperture 
or at one of the comers. In the present research, 
the origin is chosen to be at the bottom left hand 
comer of the camera aperture. An image is de- 
fined as an irradiance distribution B(x p , Zp ) in the 


_ f{x - Xp) 

p (y-yo) 


( 1 ) 


_ /(* ~ *o) 

p (y - yo ) 


( 2 ) 


Here, / is the imaging device focal length. As in 
the research of Horn and Schunck 3,4,8 , the follow- 
ing assumptions are next invoked : 

1. The underlying image irradiance distribution 
is smooth in the sense that first spatial par- 
tial derivatives exist everywhere on the image, 
except at finite number of regions. 
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2. The perceived change in image irradiance at 
each point in the image plane is entirely due 
to the motion of the imaging device, and not 
due to changes in the scene reflectance. 


These assumptions lead to the expression 



(3) 


It is to be emphasized that higher order deriva- 
tives of the irradiance may also be equated to zero. 
This fact may be exploited to resolve navigational 
anomalies in certain situations. However, the re- 
sulting navigational scheme will be more complex. 
This issue will not be pursued any further in this 
paper. Next, expanding equation (3) yields 


dE dE . BE . 
ot + dx p Xp + d7 p Zp - 0 


(4) 


Expression (4) is the Horn-Schunk optical flow 
constraint equation. The apparent velocity of the 
image x p ,z p in this expression are called optical 
flow components. 

Next, differentiating the perspective projec- 
tion equations (1), (2) with respect to time and 
noting that the objects sire fixed relative to the 
earth fixed frame leads to the equations : 


. _ [Zpilo - fig] 

p (v-vo) 


(5) 


. _ [*pV o ~ /io] 

p (y - yo) 


( 6 ) 


Substituting these relations in the optical flow con- 
straint equation results in 


that anomalies of various kinds can arise while ap- 
plying this equation to real image sequences. For 
instance, points in the image at which the spatial 
irradiance gradients are discontinuous, the expres- 
sion (7) may not have a solution. It is assumed 
here that such points can be edited-out of the im- 
age sequence using additional criteria. 

Before proceeding to use equation (7), the dis- 
crete nature of the imaging process needs to be 
examined to determine the acceptable spatial and 
temporal data rates required for obtaining a re- 
alistic estimate of the scene depth. Clearly, the 
trade-off between spatial and temporal sampling 
intervals must be examined before synthesizing a 
numerical algorithm. The simplest way to carry- 
out this analysis is to replace the partial deriva- 
tives in expression (7) by finite differences. This 
yields 


AEi 

At 


' 2, . ,. , 


+ ~ {z P y 0 - fz 0) 


Az p 


J {y - yo) 


( 8 ) 


Here A E\ is the difference in the irradiance at a 
point x p ,Zp at the begining and end of temporal 
sampling, AE 2 is the irradiance difference varying 
x p while holding the time t and z p fixed, and A E 3 


is the irradiance difference varying z p while hold- 
ing t and x p fixed. At, Ax p , A z p are the spatial 
and temporal sampling intervals. To ensure a nu- 
merically well-conditioned system, the spatial and 
temporal irradiance differences should be more or 
less equal. Imposing this requirement in (8) leads 
to : 


dE 

3t 


\dE • /• x 

Q^~( X PV0 ~ J X 0 ) 


dE, . 

+ ^(W - fz 0) 


(y - yo) 


(7) 


Expression (7) is the optical navigation equation. 
If the spatial and temporal partial derivatives of 
the image irradiance at any instant are known, 
equation (7) may be used to solve for the depth 
(y ~ Vo) since the vehicle velocity components 
* 0 , yo, io and a particular pixel location in the im- 
age plane x p , z p are known. It is important to note 


At = Ax p Az p (y - y 0 )/N r (9) 

where 


N r = ~ 


&z P (x p yo - fio) 


+Ax p (z p y 0 - fz 0 ) 


( 10 ) 


Thus, if the required depth sensitivity, vehicle ve- 
locity components, camera focal length, pixel lo- 
cation and the spacing between pixels along the 
Xp — Z p directions are known, desirable temporal 
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sampling interval can be computed using this ex- 
pression. Since the sampling interval given by ( 9 ) 
depends on the pixel location, different sampling 
intervals will be specified at various parts of the 
image. For instance, if the vehicle were travelling 
along the Y direction in the earth fixed frame, the 
sampling interval at the center of the image can 
be very large, while the sampling interval at the 
image periphery should be small. In this situa- 
tion, the smelliest sampling interval may be used 
for the entire image to simplify calculations. It 
is important to remember that the expression ( 9 ) 
provides only an upper bound on the sampling in- 
terval. Normally, a sampling interval between one 
fifth to one tenth this value should be employed. 

Once the spatial and temporal sampling in- 
tervals are selected, the only remaining issue is the 
method for computing partial derivatives. There 
are several techniques available for computing par- 
tial derivatives in the literature 13 . Most of these 
techniques suggest data smoothing prior to deriva- 
tive computations. This is an essential step be- 
cause derivative computation is an inherently noise 
amplifying process. In the following, a technique 
based on optimal control theory will be advanced 
for the computation of spatial and temporal par- 
tial derivatives. 


scheme is to produce a relatively noise free esti- 
mate of spatial and temporal partial derivatives. 
To meet this objective, define the instantaneous 
error between the model and the actual image as 

r(ip, Zp , f) = E(x p , z p , f ) Em {xp, Zp , f) (11) 

It is desired to compute the partial derivatives of 
the image model dE M /dt, dE M /dx p , dE M ldz p us- 
ing the measured irradiance and the error between 
the image model and the actual image. The im- 
age model will be continuously updated to closely 
track the actual image. Then the partial deriva- 
tives can be computed as 


6Em 

dE 

dr 

( 12 ) 

dt 

: dt 

di 

6Em 

dE 

dr 

(13) 

dip 

dxp 

dXp 

QEm 

dE 

dr 

(14) 

II 

: ^ 
<0 

dzp 

dZp 


In order to develop estimates of the error par- 
tial derivatives, consider the following optimiza- 
tion problem 


3 Estimation of Partial 
Derivatives 


min r rrwj 

r ( x p,Zp,t) J to J XpO J Zp 0 l 


dr 

[dx 


12 


P*l 



The success of the proposed navigation depends on 
estimating a consistent set of image irradiance par- 
tial derivatives. The finite difference technique has 
been the main approach employed by previous in- 
vestigators while using the optical flow constraint 
equation. While the simplicity of this approach 
is attractive, it has no noise attenuating proper- 
ties. This is a serious limitation while dealing with 
real image sequences. An alternate approach for 
partial derivative estimation based on the integral 
least square error criterion will be developed in the 
following. 

Let the model of the image under considera- 
tion be Em{zpi z p*t) and the actual image is given 
by E(x P i z pt t). The objective of the estimation 


(15) 

with given constants a, £, 7 . These constants can 
be chosen to regulate the relative magnitudes of 
the error partial derivatives. A large value would 
tend to decrease the partied derivatives implying a 
nearly uniform error distribution, while a smaller 
weight would permit larger error variations in the 
image plane. Additioneilly, these weights can be 
used to correct for errors in the camera optics by 
constraining the partial derivative magnitudes in 
certain directions. In the integral (15), t 0 is the 
initial time, x p q and x p f define the left and right 
edges of the image, and z p 0 and z pf are the bottom 
and top edges of the image. 
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The problem defined in (15) is a least square 
error estimation problem involving three indepen- 
dent variables. The necessary condition for opti- 
mality is given by 9 


r 


d 2 r Q 2 r d 2 r 

a d^ p ~ /3 d^~ 1 di 2 


= o 


(16) 


Equation (16) is a second order linear partial dif- 
ferential equation and can be solved for certain set 
of boundary conditions 11 . In the present situation, 
the image is rectangular and the error is assumed 
to decrease from the left edge to the right edge 
and from bottom edge to the top edge. It is desir- 
able that the tracking error decreases monotoni- 
cally over time. Additionally, the error is required 
to approach zero for sufficiently large t, x p , Zp. 

An approach for obtaining the solution to (16) 
is to assume that the image error has a decompo- 
sition of the form 


r = r 1 (i)r 2 (ip)r 3 (2p) (17) 

The only justification that can be advanced for 
using such a decomposition is that the resulting 
system of equations can be solved in closed form. 
Substituting this decomposition in the necessary 
condition for optimality (16) permits one to con- 
vert the second order partial differential equation 
into three second order linear ordinary differential 
equations of the form 


d?ri 

C 1 

(18) 

a 2 

= — ri 

7 

d 2 r 2 

c 2 

(19) 

dz 2 

= — r 2 
a 

cPr 3 

C 3 

(20) 

dz 2 

= jr 3 

Here, Ci^C^Cz are three constants 
the algebraic equation 

which satisfy 


1 - C 1 - C 2 - C z = 0 (20) 

Monotonically decreasing error in both spatial and 
temporal directions requires 

Ci > 0, C 2 > 0, C 3 > 0 (21) 


In this case, a closed form solution for equation(16) 
can be shown to be 


r — J{ e -l^(^o)^b(x p ~xpo)-\-c(z p -z p o)] 


( 22 ) 


where 


a = yfc^h, b = y/cJZ, c = y/cI/0 (23) 

and AT is an arbitrary constant. In expression (23), 
any two of the three constants C\ y C*i,Cz can be 
arbitrarily chosen, the third one being determined 
by the constraint equation (20). 

Differentiating expression (22), the spatial 
and temporal partial derivatives of the error can 
be computed as 


dr 

dx p 


= —br 


dr 

dz p 


= —cr 


(24) 

(25) 


dr 

di 

Thus, at any pixel in the image plane, the spatial 
and temporal partied derivatives of the error can 
be computed if the spatial find temporal error be- 
tween the image model and the actual image are 
known. Next, recognizing the fact that the image 
is spatially discretized, the image model may be 
updated at each pixel by integrating its temporal 
partial derivative with appropriate initial condi- 
tions. For instance, if an Eider integration scheme 
were employed for the temporal partial derivative, 



Em{U+i) = Em(U) + h 


l ^+a(E-E M ) 


t=ti 

(27) 


In this expression, h is the integration step size, 
normally chosen to be the image frame rate. Ex- 
pression (27) is the image irradiance tracking equa- 
tion which ensures that the image model is close 
to the actual image. The actual image partial 
derivatives dE/dt,dE/dx p ,dE/dz p are computed 
using finite difference approximations. In image 
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sequences with relatively unchanging spatial irra- 
diance gradients, these may be assumed zero, let- 
ting the tracking process to compensate for this 
value. 

The foregoing analysis may next be used to 
synthesize a navigation scheme. At a particular 
pixel in the image, the spatial and temporal par- 
tial derivatives of the image can be computed us- 
ing expressions (24)-(26). The temporal error is 
then integrated with the initial condition on the 
image model irradiance to obtain the image model 
for the next time step. The temporal sampling 
interval for this image irradiance tracking scheme 
is chosen based on equation (9) and the values of 
Ci, 7. Various constants involved in the compu- 
tations can be selected based on the quality of the 
tracked image. A schematic block diagram of pro- 
posed the navigation scheme is given in Figure 2. 

4 Results and Discussions 

The navigation algorithm was implemented in a 
simulation employing a piece-wise planar scene. 
The piece- wise planar scene is made up of joined 
rectangular planes of equal height with a contin- 
uous irradiance distribution across the surface. A 
horizontal section through this simulated scene is 
given in Figure 3. Corresponding irradiance dis- 
tribution is given in Figure 4. 

Figures 5 shows the imaging device trajectory, 
and both the actual and estimated scene depths 
at a pixel using the present navigation scheme. 
The anomalies that arise at image irradiance gra- 
dient discontinuities may be observed in this fig- 
ure. These anomalies arise due to the fact that 
the image model lags behind the actual image. As 
a result, at points where the irradiance gradients 
are discontinuous, the spatial and temporal deriva- 
tives do not change simultaneously. While this lag 
is not too serious in regions of small irradiance 
change, it becomes increasingly inaccurate when- 
ever the irradiance gradients undergo jumps. 

To a degree, these anomalies can be corrected 
by decreasing the temporal sampling interval. Al- 
ternately, one may either use a higher-order irra- 
diance tracking scheme or process the computed 
range estimates through a low-pass filter. The lat- 


ter approach was employed in the present study. 
Figures 6 and 7 show the depth estimates obtained 
using range filters with two different time con- 
stants. This investigation reveals that the pro- 
posed navigation scheme can be made to work sat- 
isfactorily in simulated image sequences. Detailed 
computations using real image sequences are cur- 
rently under way. 

5 Conclusions 

An image-based passive navigation scheme for ro- 
torcraft was described. This scheme estimates the 
scene depth by tracking the image irradiance us- 
ing a proportional feedback law. Feasibility of this 
approach was explored using a simulated image se- 
quence. 

The contributions of the present research are 
the following: 

• Derivation of the image based navigational 
equation combining the motion parameters 
with the spatial and temporal irradiance gra- 
dients. 

• Formulation of the range estimation problem 
as a distributed parameter observer problem 
with spatio-temporal irradiance distribution 
as an input. 

• Feedback solution of this problem using the 
Calculus of Variations. 

This scheme is currently being evaluated us- 
ing real image sequences. Work on a navigation 
scheme using higher-order tracking schemes and 
multiple image sequences is also underway. 
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ACTUAL MAGE SEQUENCE 



Flfl. 2. Image-Based Passive Navigation Scheme 
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3. Ranging Based on Taylor Series Expansion 


The ranging scheme described in the previous chapter suffers from a major disadvantage, viz., 
since the optical flow constraint equation relates temporal irradiance partial derivative with the spatial 
irradiance partial derivatives and optical flow, it is not clear how it can be modified for stereo image 
sequences. Furthermore, since the images are discrete, the concept of a temporal partial derivative is 
artificial. 

To alleviate these conceptual difficulties, multi-dimensional Taylor series approximation of the 
correspondence hypothesis is next introduced. Given an image E(x p , y p ), define the gradient vector 
g and the Hessian matrix H as: 


’ dE/dx p ' 

, H = 

d 2 E/dx p 2 

d 2 E/dx p dy p 

dE/dy-p 

d 2 E/dy p dx p 

d 2 E/dy p 2 


For a pair of images E \ , Ei, the correspondence hypothesis may be stated as: 


(3.1) 


E\(x v , j/p) = E 2 (x p + Ax p2 , y p + Aj/ p2 ) 


(3-2) 


E 2 (x p , j/p) — E\[x p -|- Axpi, y p T Aj/pi) 
Defining the disparity vector r, 


(3.3) 


Axp 

A j/p 


the correspondence hypothesis can be expanded in a Taylor series as : 


(3-4) 


£] - E 2 = g r 2 r 2 + ^r T 2 H 2 r 2 + (3.5) 

E 2 - E x = g T jr, + + (3.6) 

Such an expansion does not require any explicit reference to time, enabling the development of an 

algorithm that is useful for both stereo and motion image sequences. Additionally, the algorithm 

performance can be improved by including higher-order terms in the Taylor series approximation. A 
benefit of this formulation is that it no longer needs the concept of optical flow, a concept that is often 
stated to be artificial in real imaging devices [4]. 

Next, the incremental perspective projection relationships at a pixel (x p , y p ) can be found to be 


\ (XpA--o - /Ax 0 ) A {y p Az 0 - fAy 0 ) /0 

Ax p = — : — : — . a 3 /p = — — - — - 7 — (3.7) 

These incremental perspective projection relationships can be combined with either (3.5) or (3.6) to 
yield depth. An alternate procedure is to assume that an object at (x p ,y p ) in the first image is the 
same as that observed at ( x p ,y p ) in the second image. This assumption leads to ri = — r 2 . In this case, 
adding and subtracting the two expansions and substituting for the incremental perspective projection 
relationships yield the ranging and range error equations as 
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e>-e 2 = - 


{ dx' r + dxl )[x ' A *° /Alo!+{ lf + ^ }ti,I ’ A " 0 “ /A! '» 1 ]^^A^ + (38) 


,dE 1 dE, 


0 = [ ( §f “ If ■ fAxo] + { w r ~ ^ ,fcA *° ■ /Aw,) . 


1 

z - z 0 - Az 0 


+ 


(3.9) 


The first equation can be used for computing the range ( z-z 0 ) while the second equation is useful 
for determining the error involved in the range calculation. If the computational resources permit, it 
may be beneficial to use higher-order terms in this series. Note that these expressions have not been 
previously obtained in the literature. 

The spatial partial derivatives necessary for the computation of range can be obtained using finite 
difference approximations. Further details on algorithm derivation and implementation is given in the 
following paper, which was presented at 1990 AIAA Guidance, Navigation, and Control Conference. 
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Abstract 

Flight vehicles operating at low altitudes such 
as rotorcraft require range determination schemes 
for locating the terrain and obstacles. The de- 
velopment of a ranging scheme combining image 
sequences from vehicle mounted passive imaging 
sensors and the vehicle motion variables obtained 
from an on-board inertial navigation system is 
described. This approach can handle sequences 
from more than one imaging device. Derivation of 
the numerical algorithm and the performance re- 
sults using a laboratory image sequence are given. 
Other applications of the proposed approach in- 
clude ranging schemes for autonomous planetary 
rovers and telerobots. 

Introduction 

The development of position determination 
schemes such as inertial navigation systems (INS) 
and global positioning systems (GPS) has now 
reached a point where flight vehicles could be lo- 
cated within a few meters with respect to pre- 
defined coordinate frames at any point on earth. 
While the accuracy of these systems are impres- 
sive, they cannot provide acceptable solutions to 
the problem of relative position determination in a 
changing or -uncertain environment. This is due to 
the fact that it is neither possible nor is it desirable 
to construct apriori data bases incorporating every 
detail of the operational environment. Use of con- 

*Member AIAA, School of Aerospace Engineering, Geor- 
gia Institute of technology, Atlanta. Mailing Address : FSN 
Branch, M.S. 210-9, NASA Ames Research Center 

Research Scientist, Member AIAA, FSN Branch, MS 
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ventional imaging Radar or the more recent scan- 
ning Laser range finders are obvious solutions to 
this problem. However, these devices require the 
radiation of electromagnetic energy which may not 
be desirable if the on-board power availability is 
low or if covert missions are contemplated. These 
factors have motivated the recent research on rang- 
ing schemes using passive imaging devices such as 
conventional or low light level television and infra- 
red imaging sensors. Indeed, passive ranging has 
emerged as a central issue in nap-of-the-earth he- 
licopter flight guidance 1,2 and autonomous plane- 
tary rover mission 3,4 * 5 . The research reported in 
this paper is primarily motivated by the current 
interest in automating the nap-of-the-earth heli- 
copter flight. 

In the past, the research activity in the passive 
ranging problem was largely driven by the robotics 
discipline 6 . Several approaches to the image based 
range determination problem have been suggested 
in recent literature 6-10 . Specifically, Reference 7 
outlines the development of an image-based recur- 
sive range determination scheme for nap-of-the- 
earth helicopter flight. In that approach, various 
features of interest in an image sequence such as 
the edges or regions of high contrast were used 
to determine the range to various objects within 
the field- of- view. Techniques such as this may 
be termed as feature-based ranging to distinguish 
them from techniques that do not explicitly use 
any features for range determination. Image-based 
ranging algorithms that do not explicitly employ 
features are termed field-based approaches . 

The number of simultaneous imaging de- 
vices employed defines another category of ranging 
schemes. For example, the task of range estima- 
tion using a temporal sequence of images from a 
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single imaging device is called the motion ranging 
problem or the Cyclopian ranging problem after 
the one eyed monster from Greek mythology. On 
the other hand, schemes employing images from 
two or more cameras simultaneously are called 
Stereo ranging schemes. 

In the robotics literature, ranging tasks ap- 
pear to include the simultaneous determination 
of vehicle motion parameters also 6,9 " 12 . In the 
case of a single imaging device, this process can 
be shown to produce an under-determined system 
and the solution can be obtained only by impos- 
ing an additional constraint. In the approach de- 
veloped in the present work, the vehicle motion 
parameters and the relative location of the imag- 
ing devices are assumed to be known, perhaps 
from the on-board inertial navigation system. The 
ranging algorithm uses this data together with the 
given image sequence to construct the distance to 
various points within the field- of- view. 

Images obtained from sensors mounted on a 
moving vehicle will exhibit irradiance changes at 
each pixel. These changes depend on the relative 
location of various objects, vehicle motion param- 
eters, relative location of the imaging devices, the 
scene surface reflectance and the location of illu- 
mination, If the surface reflectance and the illumi- 
nation are assumed to remain constant during the 
knaging process, then the observed image irradi- 
ance changes at each pixel are entirely due to the 
relative location of various objects within the field- 
of-view and vehicle motion parameters. Addition- 
ally, if all the objects within the field of view are as- 
sumed fixed with respect to an inertial frame, then 
the irradiance change at each pixel location is in- 
dicative of the range to these objects. In this case, 
if the vehicle motion parameters and the imaging 
device constants are known, it is possible to deter- 
mine the relative location of various points within 
the field of view. The ensuing analysis will deal 
with monochromatic images. The data available 
in such an image is a distribution of irradiance 
specified as a function of imaging device aperture 
coordinates. Cartesian image coordinate system 
will be used here since most imaging devices used 
currently are rectangular. 

All the ranging algorithms reported in the lit- 


erature have their basis in the so called correspon- 
dence hypothesis . According to this hypothesis, 
if it were possible to establish the correspondence 
between various objects in a pair of images and to 
measure the displacement of these objects on the 
image plane, then the perspective projection ge- 
ometry can be employed for computing the range 
to these objects. The object displacement in the 
image plane is sometimes termed as disparity in 
the machine vision literature. In feature-based 
ranging techniques, the disparity is computed only 
between a few objects of interest in the image 
plane. On the other hand, in field based tech- 
niques, the disparity is determined at every point 
in the image plane. 

In the field-based scheme discussed in this 
paper, the correspondence hypothesis is approx- 
imated using multi- dimensional Taylor series. In 
the special case where a temporal image sequence 
is considered, the lowest degree correspondence 
equation turns out to be the optical how constraint 
equation 6 * 13 * 14 . In Reference 8, with the assump- 
tion that the objects within the field-of-view were 
fixed, this equation was combined with perspec- 
tive projection geometry and the inertial naviga- 
tion system data to derive a direct ranging equa- 
tion. This ranging equation related the image irra- 
diance spatial partial derivatives with the distance 
to various objects within the field-of-view. An ex- 
pression for determining the spatial sampling in- 
tervals required for satisfactory performance of the 
proposed scheme was also derived. The central re- 
quirement in using this approach is the availabil- 
ity of spatial and temporal partial derivatives of 
the image irradiance. An optimal approach for 
estimating various partial derivatives was synthe- 
sized using the Calculus of Variations 15 . The per- 
formance of the ranging scheme was then demon- 
strated using a simulated image sequence. 

A field-based range determination technique 
was outlined in Reference 8. The chief accomplish- 
ment there was the derivation of an Optical Rang - 
ing Equation . This paper presents an improved 
version of the ranging equation together with a 
Range Error Equation . These expressions are use- 
ful for cyclopian image sequences as well as stereo 
image sequences, and require only the spatial par- 
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tial derivatives of the image irradiance. The per- 
formance of this ranging scheme is demonstrated 
using a laboratory generated image sequence. Al- 
though the algorithm is capable handling the vehi- 
cle rotation and translation, the discussions in this 
paper will be limited to the translational motion 
case. 

The Ranging Scheme 

To avoid the complications that can arise 
while considering an inverted image as is encoun- 
tered in real imaging devices, the frontal image 
plane 16 representation will be used for image de- 
scription. Consider an image coordinate system as 
shown in Figure 1. In this figure, the major axis 
of the image plane is designated as the X p axis 
and the minor axis is labelled the Y p axis. The 
Z p axis is normal to the image plane, aligned with 
the axis of the lens. The origin of this image plane 
coordinate system is located at the center of the 
aperture. In this coordinate system, an image may 
be defined as an irradiance distribution E(x p , y p ) 
in the X p , Y p plane, with E being the irradiance 
specified on a gray scale at a point x py y p on the 
image plane. It is assumed in the ensuing analy- 
sis that the image irradiance distribution is con- 
tinuous everywhere in the image plane except at 
finite number of regions. Thus, for the purposes 
of the present analysis, effects such as spatial im- 
age discretization and gray scale quantization are 
ignored. 

The lens center is assumed to be located at 
the origin of the vehicle body axis system and the 
lens axis is aligned along this coordinate system. 
The vehicle body axis system may be related to 
an inertial frame through yaw, pitch, roll rota- 
tions <f> and three translations £o, yo, ^o- For 
the present research, the vehicle attitudes are as- 
sumed to be zero find the on-board inertial navi- 
gation system is assumed to produce estimates of 
Let se,y , z be the location of an object 
within the field-of-view, specified with respect to 
the inertial frame. The coordinates corresponding 
to this object on the image plane is x p , y p . 

The point x py y p corresponding to the object 
in the image plane is related to its location with 


respect to the inertial frame x,y y z through vari- 
ous translations and rotations of the vehicle body 
frame, together with the perspective projection. 
Effectively, the perspective projection collapses a 
three dimensional scene into a two dimensional 
plane. The objective of any ranging scheme is to 
invert this projection. The development that en- 
ables such an inverse transformation is the corre- 
spondence hypothesis . 

The correspondence Hypothesis 

The fact which allows the recovery of three dimen- 
sional geometry from a two dimensional scene is 
that if one has at least two views of a scene ob- 
tained from different vantage points, then the dif- 
ference between the location of the corresponding 
objects in the two image planes is a measure of 
the range. The apparent object displacement in 
these two images is called disparity . The farther 
the objects are with respect to the imaging de- 
vice, the lesser would be the observed disparity. 
For instance, objects that are very far away would 
appear at nearly the same location in the two im- 
ages, while close objects will exhibit a large dispar- 
ity. If the disparity of the corresponding objects in 
an image sequence could be measured, then equa- 
tions for perspective projection could be used to 
compute the scene depth. 

Let Ei(x py y p ) be the image of a scene and 
Ei{x py y p ) be another view of the same scene from 
a different vantage point. Assuming that the two 
images contain all the objects of interest and that 
to a large degree, the irradiance of the scene is 
independent of the vantage points from which the 
images were obtained, one may write the relations 

E\{x py y p } = E 2 (x p + Ax P 2 , y p + Ay p2 ) (1) 


E 2 (x p ,y p ) = Ei{x p + Ax p i,y p + Ay pl ) (2) 

The equations (1), (2) state that an object at 
the point a? p , y p in the first image will appear at 
the point x p + A z p i,y p + Ay p i in the second im- 
age. Similarly, an object located at the x p ,y p 
point in the second image will appear at the point 
xp + Aa; P 2 , y p + A y p2 in the first image. Here, 
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Aa; p i, Ay pX , Az p 2 , Ay P 2 are the disparities between 
the corresponding objects in the two images. The 
depth information is encoded in these disparities. 
It will be demonstrated subsequently that if these 
quantities can be determined, then the location of 
the objects with respect to the inertial frame can 
be computed using perspective projection geome- 
try. 

Note that the correspondence hypothesis is 
not valid at every point on the two images. De- 
pending upon the relative locations of the various 
objects within the field of view and the relative lo- 
cations of the imaging devices, the correspondence 
hypothesis will hold good only in certain regions 
of the image. This is due to the various occlusions 
that can arise from the angle of viewing and the 
camera aperture dimensions. Note that if more 
than two images of the same scene were available, 
then it is possible to write correspondence equa- 
tions such as (1) and (2) between every pair of 
them. 

The correpondence hypothesis may be satis- 
fied in various ways. In general, the task of satis- 
fying the correspondence leads to an n x m search 
for 2 X n x m quantities, with n being the number 
of pixels per column of the image and m being the 
number of pixels per row of the image. That is, 
at each pixel on the image, one searches for the 
disparities along the two directions between two 
images. Although additional information is avail- 
able to limit the search domain, this process is 
impractical and computationally ill conditioned. 

A family of techniques called the feature 
based methods 7 attempt to decrease the complex- 
ity involved in satisfying the correspondence hy- 
pothesis by first extracting the features of inter- 
est in the image such as edges or regions of high 
contrast and then establishing correspondence be- 
tween them. This step reduces the correspondence 
task by a significant amount. Pixel by pixel cor- 
relation has also been employed for satisfying the 
correspondence hypothesis 17 . 

In the present paper, an alternate approxi- 
mation for satisfying the correspondence hypothe- 
sis will be advanced. In this technique, the corre- 
spondence equations (1) and (2) are first expanded 
using a Taylors series, to yield 


Ei(x p ,y p ) - E 2 (z p ,y p ) = — 2 ^” yp) A Zp2 


dx v 

dE 2 (x p ,y p ) 

+ % Ay ’ 1 + 


( 3 ) 


E 2 (x P ,y P ) - E^x^yp) 


d-Ei(s P ,y P ) 

dx p 


Az p i 


+ dy r Ai " 1 + 


( 4 ) 


Expressions (3) and (4) relate the irradiance 
spatial partial derivatives to the irradiance change 
and the image disparities. The second order terms 
in these expansions will typically involve the Hes- 
sian matrix, while the third term can only be ex- 
pressed as a tensor. Truncating this series pro- 
duces various orders of approximation. Clearly, 
the linear approximation is preferable since it leads 
to a simple algorithm. Hence, in all that follows, 
the nonlinear terms in the Taylor series expansion 
will be dropped. 

In the special case where the disparities 
Azp, A j/p arise from vehicle motion yielding appar- 


ent image velocities u, v at a pixel in the image, i.e. 
Ax p = tiA t, Ay p = vAt, the equations (3) and (4) 
with higher order terms dropped is the optical flow 
constraint equation given by Horn and S chunk 13 . 
In this case, the velocity components u, t; are called 
the optical flow components. In the present pa- 
per, this characterization is avoided because the 


disparities Ax p \, Ay p i, Ax p2t Ay p2 may arise from 
sources other than vehicle motion. Moreover, even 
in such temporal image sequences, it has been sug- 
gested that the optical flow concept is artificial 13 . 
By adopting the point of view proposed in the fore- 
going, equations (3) and (4) may be used not only 
for temporal image sequences such as those en- 
countered in cyclopian vision problems but also for 
stereo image sequences or images obtained from 
multiple cameras. 

In the next section, it will be demonstrated 
that these expressions can be used together with 
the geometry of perspective projection to obtain 
a ranging equation. That analysis will also lead 
to the derivation of an equation for computing the 
range error. 
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The Optical Ranging Equation 

Assuming that the lens center is located at the 
origin of the vehicle body coordinate system and 
aligned with it, the relation between an image 
point Xp, y p and the position of the correspond- 
ing object with respect to the inertial frame x, y, z 
is given by the perspective projection equations 

X P = f( x - x o)/{z - *o) (5) 

y P = f{y-y o)/(*- z 0 ) (6) 

Here, / is the imaging device focal length, and 
zoiVoiZo are the location of the vehicle with re- 
spect to the inertial frame. Note that if the cam- 
eras are offset from the origin of the vehicle body 
axis system and not aligned with it, additional co- 
ordinate transformations will have to be incorpo- 
rated in the analysis. 

Next consider a translation of the imaging de- 
vice by a small distance Ax 0 , Ay 0 , Az 0 in the pos- 
itive direction. These variables could also be in- 
terpreted as the displacement of the second imag- 
ing device relative to the first as in stereo rang- 
ing schemes. It is assumed that the location of 
the object point s,y,z remains unchanged during 
this motion. In this case, the object located at the 
point s, y, z in the scene, and observed at the point 
x p , y p in the previous image plane would now ap- 
pear at the point x p + Ax p , y p + Ay p on the image 
plane, i. e 


x p + Ax p = f(x - s 0 - As 0 )/(z - zq- Az 0 ) (7) 


Vp + Ay p = /(y - y 0 - Ay 0 )/(z - zq - Az 0 ) (8) 

Subtracting (5), (6) from the expressions (7), 
(8), results in 


a- _ (x p Azo- fAx 0 ) 

LXX p — “ 

z - z 0 - Azo 


( 9 ) 


a „ _ (vp Az o - f A v o) 

i\y p — — 

z - zq- Az 0 

Note that if the direction of motion is changed, the 

resulting displacement in the image plane can be 



computed by changing the signs of the variables 
Ax 0 ,Ay 0 ,Azo' 

Expressions (9) and (10) may next be used to 
eliminate the variables Ax pl , Ay pl , Ax p2 , Ay p2 in 
the correspondence equation pair (3), (4). At each 
pixel location x p , y p this yields : 


f dE 

Ei - E 2 = j ^-^[zpAzo - fAx 0 ] 

+ ^ A ^ /Aso] }r -z!-A^ < u > 

f dE 

E2 - Ei = -|-^^ : [xpAzo - fAx 0 ] 

+ *jr 1 [y P Az o-f A yo}} 1 . ( 12 ) 

oy p J z - z 0 - Azo 7 

Adding and subtracting the expressions (11) and 
(12), one may rewrite these two equations as 
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f 9 Ei dE 2 , r 

{ dZ + - /AxcI 


+{ ^ + i> A2 »- /Aml 


z - z 0 - Az 0 
(13) 




dE 2 

dx p 


}[s p Az 0 - fAx 0 ] 



qe 2 

dy P 


}[ 2 / P Azo-/Ay 0 ] 


1 

Z — Zq — Az 0 


(14) 

The expression (13) is the Optical Ranging equa- 
tion relating the irradiance difference at the same 
pixel locations in the two images with the spa- 
tial partial derivatives in the image planes and the 
scene depth. This expression has not been previ- 
ously reported. 

The expression (14), on the other hand, pre- 
dicts the degree of error in approximating the cor- 
respondence equations using the first two terms of 
the Taylor series expansion. This is a consequence 
of the fact that if the image irradiance distribution 
was indeed linear, the spatial partifil derivatives 


17 



will remain the same in both images E\ and Ei. 
Interestingly, the expression (14) also reveals that 
for the same contrast level, image-based range es- 
timation is more accurate for objects that are near. 
Assuming that the irradiance differences ( E\ — E 2 ) 
are nonzero, the expression (13) can be manipu- 
lated to the form 


* dE 9 E 

z = z 0 + Azo+ + -^-}[x p Azq - f Axq) 


+{ 


dE x dE 


+ 


dy P dy 


-}[y P Az 0 - fAy 0 ] 


,J J 2 (E 1 - E 2 ) 

(15) 

Expression (15) can be used to compute the scene 
depth z if the image irradiances E l5 E 2 , camera 
parameters /, s p , y p , Ax 0s Aj/O) Azo and the spatial 
partial derivatives 9Ei/dz p , dE x jdy p , dE 2 /dz p , 
dE 2 /dy p were available. 

Next, the se,y coordinates can be determined 
using the inverse of the perspective projection 
equations together with the optical ranging equa- 
tion as follows. 


A zq dE\ 9E 2 

x - x ° + ~1~ + aT, * aZ' }[tpA2 ° " /Ax ») 


+< H + f | }[! " A2 »- /A "-1 


2f(£, - B,) 

(16) 


y p Azo \ f dEi dE 2 .. 
y - »o + -±j- + [{j- + j- }{zpAz 0 - fAxo] 


f dE i 
+{^~ i + 


dy p dy. 


9E2 }[y P Azo-fAyo] 


Vp 


I 2 f{E x - E 2 ) 

(17) 

In the preceding development, it is important 
to note that anomalies of various kinds can arise 
while applying the ranging equation to real image 
sequences. This is because of the fact that this 
development assumes that the first partied deriva- 
tives exist and that at least locally, a linear ap- 
proximation is valid for the image irradiance. Ad- 
ditionally, in regions of near uniform irradiance, or 
in regions where the spatial irradiance gradients 
are discontinuous, the expressions (15)-(17) may 


not have a solution. It is assumed here that such 
points can be edited-out of the image sequence by 
the use of an appropriate threshold. 

The Computational Algorithm 

The first step in the implementation of the rang- 
ing algorithm is the computation of image spa- 
tial partial derivatives. Several techniques are 
available for computing partial derivatives in the 
literature 6 . A prominent one amongst them is the 
Fourier transform technique. Although fast algo- 
rithms are available for their computation, the cal- 
culations are far from routine. In Reference 8, a 
partial derivative estimation scheme based on opti- 
mal control theory was advanced. Since this tech- 
nique has not been completely evaluated, a simple 
central differencing scheme will be used to com- 
pute the spatial partial derivatives in the present 
research. 

According to the central difference scheme, if 
i is the pixel location along the z p direction and j 
the pixel location along the y p direction, with the 
distance between pixels being 6x,6y , the spatial 
partial derivatives are can be approximated by the 
relations 


dE _ E(i + 1) - E{i - 1) 
dx p ~ 2 5x 


(18) 


8E + 

9, r ~ Wy < 19 > 

The partial derivative estimates given by the 
expressions (18) and (19) are rather crude, and do 
not exhibit any noise attenuating characteristics 
at all. Research on more advanced partial deriva- 
tive computation schemes are currently underway. 
Once the spatial partial derivatives are computed, 
equations (15), (16), (17) can be used to determine 
the location of various objects within the field-of- 
view. 

A block diagram of the ranging algorithm is 
shown in Figure 2. Two views of the candidate 
scene are first used to compute the spatial deriva- 
tives of the images. Next, the second image is sub- 
tracted from the first to obtain irradiance changes 
at each pixel. In the regions where the irradance 
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change is close to zero, the scene depth is not com- 
puted, since the ranging equation is not valid in 
these regions. Instead, these points are assigned 
the maximum scene depth. The spatial gradients 
and the irradiance changes are then combined with 
imaging device coordinates in the given inertial 
frame x 0 , y 0j Ax 0 , Aj/o, A z 0 to obtain the x, y, z 
coordinates. 

When a new image becomes available, the first 
image in the sequence is dropped and these cal- 
culations are repeated. In some image sequences, 
the estimated range can be inaccurate due to noise, 
digitization errors or computational errors. In this 
situation, one may collect several images of the 
scene from the same vantage point. These images 
may then be averaged to yield relatively noise-free 
range estimates. An alternate approach would be 
to apply a defocusing operator 6 to the image se- 
quence before carrying out the range calculations. 

Results and Discussions 

Several laboratory image sequences were used 
for demonstrating the performance of the ranging 
algorithm. A pair of images from this laboratory 
image sequence is shown in Figures 3 and 4. The 
scene consisted of a dark wall with a table in the 
foreground. Two computer monitors, an H shaped 
object and three pencils were then placed on this 
table. This image sequence was obtained using 
a CCD (Charge Coupled Device) camera with a 
field-of-view of about 45 degrees and a focal length 
of about 6 mm. The gray scale consisted of 256 lev- 
els corresponding to 8 bit digitization of the CCD 
image. The CCD output was digitized and pro- 
cessed on a workstation. 

Gray scale coded spatial partial derivatives 
corresponding to the images in Figure 3 find 4 are 
shown in Figures 5 and 6. Note that these are the 
average of the irTadiance partied derivatives in the 
two images. In order to produce this representa- 
tion, the gradients are first normalized such that 
they are in the range ±128. The result is then 
biased by 128 such that the level black represents 
the maximum negative gradient, while the white 
level represents the maximum positive gradient. 
Figure 7 gives the irradiance difference between 


two images shown in Figure 3 and 4. Once again, 
the differences are normalized and scaled such that 
the largest negative value would correspond to the 
black level find the largest positive value would 
correspond to the white level. Based on the im- 
age digitizer characteristics, in regions where the 
irradiance differences were smaller than 15 levels 
of gray, range computations were not carried out. 

Range calculations were next carried out us- 
ing this image irradiance data. Although a 3-D 
representation of the result was desirable, no sat- 
isfactory computer packages were available. In- 
stead, the range data was once again coded in 
terms of a gray scale. Figure 8 shows the result 
of applying the ranging equation to the present 
image sequence. In this figure, objects that are 
at zero range are represented by the black color, 
while the points that are far, or the points at which 
range data is not available are indicated in white. 
It may observed from this figure that the relative 
location of various objects in the image sequence 
are clearly recovered. The range computations ap- 
pear sparse because of the threshold employed at 
the input. Lowering this threshold will result in 
a denser range field, but will also contribute to 
several inaccurate range values. 

To a degree, these anomalies can be corrected 
by decreasing the image sampling interval. Al- 
ternately, one may either use an improved partial 
derivative estimation scheme or implement a rang- 
ing algorithm based on higher-order correspon- 
dence equation. The present investigation reveals 
that the proposed ranging scheme works satisfac- 
torily in laboratory image sequences. Work is cur- 
rently underway to apply this technique to image 
sequences obtained from a helicopter flight. 

Conclusions 

An image-based passive navigation scheme for 
rotor craft was described. This scheme estimates 
the scene depth using irradiance differences and 
spatial partial derivatives of an image sequence 
together with imaging device parameters. Perfor- 
mance of the ranging algorithm was demonstrated 
using a laboratory image sequence. 
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The contributions of the present research are 
the following: 

• Development of an image correspondence 
equation valid for both cyclopian and stereo 
ranging schemes. 

• Derivation of the Optical Ranging Equation 
combining the vehicle and imaging device pa- 
rameters with irradiance differences and the 
spatial irradiance gradients. 

• Derivation of a Range Error Equation to aid 
in evaluating the results generated by the 
ranging equation. 

• Demonstration of the performance of this al- 
gorithm using a laboratory image sequence. 

Improvements on this ranging algorithm is an 
on-going research activity. 
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Fig. 4. Second Image of the Laboratory Scene , ^*8* ®* ® ra y Scale Coded Spatial Par- 

tial Derivative of the Laboratory Scene in 
the Vertical Direction 



Fig. 5. Gray Scale Coded Spatial Par- 
tial Derivative of the Laboratory Scene in 
the Horizontal Direction 


Fig. 7. Gray Scale Coded Irradiance 
Difference Between the Two Images at Each 
Pixel 
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4. A Partial Derivative Estimation Method and its use in 

Vision-Based Ranging 


In the previous chapter it was demonstrated that vision-based ranging can be accomplished if 
the image irradiance spatial partial derivatives were known together with camera parameters. Due to 
the noise encountered in real image sequences, it is desirable to carry out noise attenuation together 
with partial derivative estimation. Techniques discussed in the literature for noise attenuated partial 
derivative estimation are all frequency domain based [2]. Such an approach is inconvenient because 
the complete image must be available before the partial derivatives can be determined. 

An alternative approach permitting derivative computation as a part of image collection was next 
pursued. The approach employs an approximation popular in the study of linear partial differential 
equations arising in mechanical vibrations. The approximation involves the factoring of a function of 
two variables as a product of two functions that vary as a function of one independent variable each. 
This factorization reduces a partial differential equation into a set of ordinary differential equation. 
Thus, the partial derivative estimation problem can be expressed as a set of total derivative estimation 
problems for which solutions already exist. Additionally, by employing causal estimators, the partial 
derivative computation can be carried out as the image data is being collected. Combining this partial 
derivative estimation scheme with the ranging equation developed in Chapter 3 then results in a fast 
method for vision-based ranging. A paper describing this research was presented at the 1991 SPIE 
Symposium on Optical Engineering and Photonics in Aerospace Sensings and is given in this chapter. 

Up to this point, the vision-based ranging formulations assumed that the camera orientation 
in the pair of images were the same. In practical situations, it is important to include both the 
translational and rotational displacements of the imaging devices. Such an extension is nontrivial 
since additional reference frames that may be translating and rotating with respect to the camera 
frame need to be included in the analysis. A paper discussing these extensions was presented at the 
1991 AIAA Guidance , Navigation, and Control Conference and is also included in this chapter. 
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A Fast Algorithm for Image-Based Ranging 
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ABSTRACT 

Image-based ranging has emerged as a critical issue in the low altitude operation of flight vehicles 
such as rotorcraft and planetary landers. These flight regimes require ranging systems for recovering the 
geometry of the terrain and obstacles for use with guidance algorithms. The development of a ranging 
equation combining image irradiance together with various order spatial partial derivatives and the vehicle 
motion parameters is discussed. The ranging equation is in the form of a polynomial in scene depth. Two- 
dimensional linear filters are then used to compute the coefficients of this polynomial to result in a fast 
image- based ranging algorithm. Performance of the algorithm is demonstrated using laboratory images. 

1 INTRODUCTION 

The problem of determining the position of various obstacles in a changing or uncertain environment 

is central to the development of guidance policies for autonomous vehicles. Imaging Radar and scanning 
•School of Aerospace Engineering, Georgia Institute of technology, Atlanta. Mailing Address : MS 210-9, NASA Ames 
Research Center 

* Member of the Professional Staff, Sterling Software 
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Laser range finders are two direct solutions to this problem. However, these devices require the radia- 
tion of electromagnetic energy which may not be desirable if the on-board power availability is low or if 
covert missions are contemplated. Recent research on ranging schemes using passive imaging devices such 
as conventional or low light level television and infra-red imaging sensors is motivated by these factors. 
Specifically, the passive ranging problem is considered a critical issue in nap-of-the-earth helicopter flight 
guidance 1,2 and autonomous planetary rover mission 3,4,5 . 

Most of the research activity in the image-based ranging problem is driven by the Robotics discipline 6 . 
Recently, several ranging approaches have been suggested in the literature 6-11 . These techniques are based 
on the intuitive notion that the images obtained from sensors mounted on a moving vehicle will exhibit 
changes in irradiance at each pixel location as a function of time. These changes depend on the relative 
location of various objects in the field-of-view, vehicle motion parameters, location of the imaging devices 
with respect to the vehicle, the scene surface reflectance and the location of illumination. If the surface 
reflectance and the illumination are assumed to remain constant during the imaging process, then the 
observed image irradiance changes at each pixel can be related to the location of various objects within 
the field-of-view and the camera- vehicle parameters. Additionally, if all the objects within the field of view 
are assumed fixed with respect to an inertia! frame, the irradiance change at each pixel location can be 
related to the range to these objects. In this case, if the vehicle motion parameters and the imaging device 
constants are known, the relative location of various points within the field of view can be determined. 
The ensuing analysis will deal with monochromatic images. A Cartesian image coordinate system will be 
used here since most imaging devices used currently are rectangular. 

In the Robotics literature, ranging tasks appear to include the simultaneous determination of vehicle 
motion parameters also 6,9,12 . If a single imaging device is used, this process produces an under-determined 
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system of equations and the solution can be obtained only by imposing an additional constraint. In the 
present work, the vehicle motion parameters and the relative location of the imaging devices are assumed 
to be known from an inertial navigation system on-board the vehicle. The ranging scheme uses this data 
together with the given image sequence to construct the distance to various points within the field-of-view. 

Image-based ranging schemes use the correspondence between various objects in an image sequence 
to determine the disparities . These disparities can be used together with perspective projection equations 
to obtain the range. In this paper, the correspondence problem is first satisfied using a multi-dimensional 
Taylor series involving disparities between the two images. Next, the disparities are eliminated in favor of 
the vehicle-camera motion parameters and scene depth using the perspective projection relationships. This 
process yields a polynomial for the scene depth at every pixel location. This equation is referred to as the 
Optical Ranging Polynomial in the present work. This equation may be thought of as the generalization of 
the optical ranging equation discussed in Reference 11. The coefficients of the optical ranging polynomial 
depend on the irradiance difference between the two images and various order irradiance spatial partial 
derivatives. A system of two-dimensional linear filters are next synthesized to compute the required partial 
derivatives. Due to the simplicity in implementing these filters, the partial derivatives can be computed 
at a high data rate, often faster than the video frame rates. The result is a fast image-based ranging 
algorithm. The performance of the ranging scheme is demonstrated using laboratory images. 

2 THE RANGING SCHEME 

To maintain consistency with the existing literature, the frontal image plane 13 representation will be 
used for image description. The various coordinate systems employed in this paper are illustrated in Figure 
1. This figure shows the image coordinate system together with the vehicle body coordinate system and 
the inertial coordinate system. The origin of the image plane coordinate system is located at the center of 
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the aperture. An image is defined as an irradiance distribution E(x p ,y p ) in this coordinate system, with 
E being the irradiance specified in a gray scale at a point x p ,y p on the image plane. The ensuing analysis 
assumes that the image irradiance distribution is continuous everywhere in the image plane except at fi- 
nite number of regions. This is a reasonable assumption because the underlying imaging process involves 
integration, a process that smoothens out the discontinuities. For the purposes of the present analysis, 
errors arising out of spatial image discretization and gray scale quantization will be ignored. 



Figure 1: The Coordinate System 


To simplify the development, the lens center is assumed to be located at the origin of the vehicle body 
axis system and the lens axis is aligned along this coordinate system. Moreover, the analysis will assume 
only pure translational motions for the camera. Although the procedure to include the rotational motion 
is rather direct, it will be developed in a more complete detail elsewhere. The on-board inertial navigation 
system is assumed to produce accurate estimates of camera-vehicle displacement x 0 ,y 0 ,Zo- In the case of 
stereo camera arrangement, it is assumed that the relative displacement between the two cameras is also 
accurately known. Let x,y,z be the location of an object within the field-of-view, specified with respect 
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to the inertial frame. The coordinates corresponding to this object on the image plane is x p ,y p . If the 
camera moves to a new position xq T Axo, yo + At/q? ~o + Azo the object point will move to a new position 
x p + Ax p , y p + A y p in the image plane. The apparent displacement in the image plane caused by the cam- 
era motion Ax p , Ay p are called the disparity. The disparity depends on the range to the object and the 
camera parameters. The objective of the ranging scheme is to use this information to determine the object 
location x,y,z. If the disparity of a corresponding object in an image pair is known, then the equations 
for perspective projection can be used to compute the scene depth. While several distinct approaches are 
available for disparity determination 6 ’ 7 , the present work will employ a multi-dimensional Taylor series 
expansion as the basis for establishing correspondence. 

Let E x {x p ,y p ) and E 2 {x p> y p ) be a pair of images of a scene obtained from two different different 
vantage points. Assuming that these two images contain all the objects of interest and that to a large 
degree, the perceived irradiance of the scene is independent of the camera location, the correspondence 
hypothesis can be expressed as 


Ei(x p ,y p ) = E 2 {x p + Ax p , y p + A y p ) (1) 

E 2 (x p , y p ) = E x (x p - A x p , y p - A y p ) (2) 

Here, Ax p , A y p are the disparities between the corresponding objects in the two images. The depth infor- 
mation is encoded in these disparities. 

Note that the correspondence hypothesis expressed in (1) and (2) are valid only in a limited sense, 
and not for all points in the two images. This is due to the various occlusions that can arise from the angle 
of viewing and the camera aperture dimensions. Note that equations such as (1), (2) can be written for 
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any number of images, provided the occluding objects are edited out. 

In the present research, the expressions (1) and (2) will be expanded in a Taylor series to enable the 
computation of the disparities. Such an expansion treats disparities as perturbations about specified pixel 
locations. The error involved in such an approximation depends on the number of terms included. Define 
the disparity vector r, the irradiance gradient vector g, and the irradiance Hessian matrix H as 


1 

a. 

H 

< 

» g = 

dE/d i p 

, H = 

d 2 E/dx p 2 

d 2 E/dx p dy p 

1 

<1 


dE/dy p 


d 2 E/dy p dx p 

d 2 E/d y 2 


Equations (1) and (2) can now be expanded using a multi-dimensional Taylor’s series to yield 

E\ - E 2 = g 2 r + ir T H 2 r + (4) 

Similarly, 

£2 - Ei = — gir+ ir T H!r - (5) 

Expressions (3) and (4) relate the irradiance spatial partial derivatives to the irradiance change and the 
image disparities. It is awkward to write the third term in this series using the vector-matrix notation 
because it contains a Tensor. In all that follows, a lower case bold face letter will denote a vector, while 
an upper case bold face letter will be used to denote a matrix. 

In the special case where the disparities Ax p , A y p arise from vehicle motion yielding apparent image 
velocities u,v at a pixel in the image, i.e. Ax p = uAt,Ay p = vAt, the equations (4) and (5) with just 
the first term on the right hand side is the optical flow constraint equation 14 . In this case, the apparent 
image velocity components u, v are called the optical flow velocities. This characterization is avoided here 
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because the disparities Ax p ,Ay p may arise from other sources also, stereo camera arrangement being an 
example. By adopting the point of view proposed in the foregoing, equations (3) and (4) may be used not 
only for temporal image sequences such as those encountered in Cyclopean vision problems, but also for 
stereo image sequences or images obtained from multiple cameras. 


Adding and subtracting the expressions (4) and (5) results in 


Ei-E 2 = - 


gl + g2 


r + -r 


H 2 - H, 


r + 


0 = 


g2 “ gl 




Hi +H 2 r + 


( 6 ) 

( 7 ) 


In the following, it will be demonstrated that these expressions can be used together with the geometry 
of perspective projection to obtain a ranging equation. The analysis will also lead to the derivation of an 
equation for computing the range error. 


The disparity vector r from equations (6) and (7) can next be eliminated in favor of other known 
quantities and the scene depth as follows. Assuming that the lens center is located at the origin of the 
vehicle body coordinate system and aligned with it, and that the camera and the vehicle are permitted 
to have only translational motions, the relationship between an image point x p , y p and the position of 
the corresponding object with respect to the inertial frame x,y,z are given by the perspective projection 
equations 

x p = f{x - x 0 )/(z - z 0 ), t/p = f(y - y 0 )/(z - z 0 ) (8) 

Here, / is the imaging device focal length, and xg, yo, zg coordinates of the vehicle location with respect to 
the inertial frame. It is assumed that the camera focal length and the components of the vehicle position 
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vector zo, j/o> *0 are known. 

Next consider a translation of the camera/vehicle by a small distance Axo,Ayo,Azo in a positive 
sense. Note that Axo, Aj/o, A^ 0 may also be interpreted as the displacement of the second imaging device 
relative to the first as in a stereo arrangement. It is assumed that the location of the object point x, y, z 
remains unchanged during this motion. In this case, the object located at the point x,y,z in the scene, 
and observed at the pixel location x p ,y p in the previous image plane would now appear at the point 
x p + Axp, y p + Ay p on the image plane. Using this fact, the disparity vector can be obtained as 

(xpA^o - / Ax 0 ) 

( y P Az 0 - fAy 0 ) 

Note that the disparity vector is linear in the inverse of the scene depth. For notational convenience, 
the inverse of the scene depth at a pixel location is denoted by the symbol 6. The quantities within the 
parenthesis in the disparity vector can be computed if the vehicle-camera displacement, pixel location and 
the camera focal length are known. Substituting for r from (9) in the Taylor series approximation (6), (7) 
yields two polynomials 


E\ — E-2 — fl] S + a? 6^ + 03 6^ -f- a 4 6 4 + ( 10 ) 

0 = 61 <5 + 62 + ^3 + 64 S 4 + ( 11 ) 


The polynomial coefficient ai depends on the vehicle motion parameters and the first spatial partial 
derivatives of the images. The coefficient <12 depends on the vehicle motion parameters and the second 
spatial partial derivatives of the two images and so on. The coefficients of the second polynomial, on the 
other hand, depends upon the differences in various order spatial partial derivatives in the two images. The 
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second equation may be thought of as a truncation error equation for the first polynomial. For example, 
if all but the first term is dropped in right hand side of equation (10), then the underlying assumption 
is that the image irradiance varies linearly in the spatial direction. In this case, the coefficient b 1 in the 
polynomial (11) is proportional to the difference between the spatial partial derivatives in the two images. 
If the two partial derivatives were exactly equal, then the error in range calculation is also zero. On the 
other hand, if they are not equal, the first term in polynomial (11) will not be zero. The magnitude of this 
term is then proportional to the error in using the first-order Taylor series approximation for satisfying the 
correspondence hypothesis. 

The solution to the polynomial (10) yields the inverse of the scene depth 6 at a particular pixel loca- 
tion. These equations have not been previously derived in the literature. In all that follows, the polynomial 
(10) will be referred to as the Optical Ranging Polynomial and the polynomial (11) will be termed as the 
Optical Banging Error Polynomial. Ihe optical ranging polynomial and the ranging error polynomial can 
be formulated at every pixel location in an image. Various numerical techniques can then be used to obtain 
6 at every pixel location in the image. This can done in closed-form for up to fourth degree approximation. 
If it is desired to use ranging polynomial of fifth degree or more, one of the several well known polynomial 
root finders can be used to obtain the solution. 

Once the inverse of the scene depth 6 at a particular pixel is computed using the optical ranging 
polynomial, the inertial coordinates of the corresponding object can be found as 

x 1 u I 1 

x — + ^~o) + i’Oi U = + 2/o, 2 = j + zq + Az 0 (12) 

It is important to observe that various kinds of anomalies can arise while using the ranging polyno- 
mial. For example, the assumption that various order partial derivatives exist everywhere in the image 
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is sometimes violated. Additionally, in regions of near uniform irradiance, or in regions where the spatial 
irradiance gradients are discontinuous, the polynomials (10)-(11) may not have a solution. To a certain 
degree, it is possible to eliminate such regions by using an appropriate threshold. 

3 The Computational Algorithm 

The first step in the implementation of the optical ranging polynomial is the computation of the 
polynomial coefficients cn,a 2 ,....a n and bi,b- 2 , ....b n . These coefficients depend on the various order spa- 
tial partial derivatives of the two images and the camera-vehicle motion parameters. In this work, the 
camera-vehicle parameters are assumed to be known from on-board inertial navigation system and geo- 
metric measurements. The remaining task is the estimation of image irradiance partial derivatives. 

Several numerical techniques are available in the literature for the computation of partial derivatives 6 . 
These range from simple techniques such as central difference schemes to very sophisticated techniques 
that adapt to reject noise and unwanted frequency components in the input signal. It is well known in the 
signal processing literature that derivative estimation is a noise amplifying process. Given an irradiance 
distribution E(x p ,y p ), the objective of the partial derivative estimators are to provide sufficiently noise 

free estimates of dE/dx p , dE/dy p ,d 2 E/dx p 2 , d 2 E/dy p 2 , d^E/dx p dy p , Most of the approaches in the 

signal processing literature formulate the partial derivative estimation problem in the frequency domain. 
This makes them unattractive in conjunction with conventional TV cameras because the image is formed 
using a sequential scanning process. 

In the following, a partial derivative estimator will be formulated directly in the spatial domain. Unlike 
the frequency domain estimators discussed in the image processing literature, the present formulation 
employs causal estimators 6 to improve speed. To a certain extent, this process compromises accuracy. 
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However, numerical experiments have shown that the error magnitude can be reduced by an appropriate 
selection of the estimator parameters. In order to simplify the development, it is assumed that the image 
irradiance distribution E(x p , y p ) can be factored into a product of two functions of one variable each. 
This allows the evaluation of two dimensional integrals as products of one dimensional integrals. In the 
following, estimated quantities will be denoted by a Thus, 

E = e\(x p ) e 2 (y p ), E = e t (x p ) e 2 (y p ) (13) 

The partial derivatives of the estimated picture E can be written in terms of total derivatives of its factors 
ci, e 2 as: 


dE_ = . de^ dE_ = . d£2 &E_ _ . m d 2 E de x de 2 

dx p 2 dx p ’ dy p 1 dy p dx p 2 62 dx p 2 ' dy p 2 €l dy p 2 ’ dx p dy p dx p dy p ““ 14 

The estimation problem can next be formulated as two, coupled, one-dimensional problems. Addi- 
tionally, it is possible to formulate the estimation problem either as a serial process or as a parallel process. 
In the former case, the processing is carried out first along the x p direction and then along the y p . In the 
latter case, the processing occurs simultaneously along both x p and y p directions. A serial implementation 
will be illustrated here since it is consistent with the image forming process. Let p be the n dimensional 
estimator state vector. The estimator dimension n is selected based on the order of partial derivatives 
to be estimated and the desired smoothness requirement. For example, if it is desired to estimate the 
first and second order partial derivatives with the requirement that the second order partial derivatives be 
smooth, n has be at least 3. Next, define a linear dynamic process evolving along the x p direction with 
the measured picture irradiance E as the input. 


d 

- — p = Ap + B£ 

ClXp 


(15) 
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The matrices A and B can be chosen to obtain the desired transmission gain and bandwidth for the x p 
spatial frequencies. A is an n x n dimensional matrix, while B is an n x 1 vector. Note that E is a scalar 
function of x p ,y p . The output equation for this linear dynamic process can be defined to yield 

r = Cp + D£,r=e 2 [e,|i 0 f (16) 

The matrices C, D should be chosen depending on the linear dynamic system (15) dimension and on the 
desired orders of the partial derivatives. 


In order to take care of the y p direction, construct a linear dynamic process in this direction with the 
first element of the vector r, viz., eie 2 as the input, i. e., 


d 

- — q = M q + N e\e 2 

ay p 


(17) 


q is an m vector, M is an m x m matrix while N is an m x 1 vector. As in the x p direction, the dimension of 
this linear dynamic process is chosen based on the desired partial derivatives and smoothness requirements. 
Once again, the matrices M and N can be chosen to obtain the desired gain and bandwidth for the y p 
spatial frequency components. The output equation of this dynamic process can be chosen as 


s = K q + Le^, s = e x [e 2 33 ^ 333 ] T 


(18) 


dy p dy p 2 

The output matrices K, L can be chosen in the same manner as the output matrices C, D in equation 
(16). With the availabilty of the vectors r,s, the estimated picture and its various partial derivatives can 
be computed. The various elements of these vectors will be denoted by lower case letters in the following. 
It can be verified that 


E — 5 i » 0 

UXr, 


E\i'2 dE _ d 2 E 
r-i ' dy p S2 ’ dx p 2 


£ir 3 (PE_ _ &E_ 

n ’ dy p 2 83 ’ dx p dy p 


£2 $2 

T 1 


(19) 
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In equation (19), it is important to ensure that rq does not become zero at any point in the image. 
This can be accomplished by adding a constant bias to the input image and subtracting this bias at the 
output. Note that such an operation does not alter the partial derivative estimates. 

It can be seen from the foregoing development that the image irradiance partial derivatives can be 
estimated as the image is being formed through the scanning process. As a result, by the time the image 
is complete, various coefficients of the ranging polynomial are available. Moreover, if a low-order ranging 
polynomial is used in the computations, the range calculations can also be completed as the image is being 
formed. 

Various steps involved in these calculations are illustrated in the block diagram given in Figure 2. 
This figure assumes that an image pair is available at the beginning of the computations. The linear 
filters provide the estimated image together with a consistent set of partial derivatives. Next, the filtered 
second image is subtracted from the first to obtain irradiance changes at each pixel. The scene depth is 
not computed in regions of near-zero irradiance change, since the ranging equation is not valid in these 
regions. Instead, these points are assigned a pre-defined maximum scene depth. The spatial gradients 
and the irradiance changes are then combined with imaging device coordinates in the given inertial frame 
zo, Ax 0 , Ayo, Azq to obtain the x,y,z coordinates corresponding to every pixel location. 


4 RESULTS AND DISCUSSIONS 

A laboratory image pair is used to demonstrate the feasibility and performance of the ranging algo- 
rithm discussed in the previous section. This laboratory image pair is shown in Figures 3. The scene 
consisted of a dark wall in the background with a table in the foreground. A soda can together with two 
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Figure 2: Vision-Based Ranging Scheme 

pencils were then arranged on this table. A wire is strung across the two pencils to explore the feasibility 
of detecting obstacles with small dimensions. These images were obtained using a CCD (Charge Coupled 
Device) camera with a field-of-view of about 45 degrees and a focal length of about 6 mm. The images 
consisted of 512 x 512 pixel matrix with the scene irradiance being discretized in 256 levels of gray. The 
stereo base line, or the lateral displacement Axq was 0.1 inches. The ranging algorithm was implemented 
on a SUN workstation in C. 


Although ranging algorithms of various order have been implemented, the results using a first degree 
ranging polynomial, together with a first-order partial derivative estimator will be presented here. The 
partial derivative estimators employed for this example are of the form 


p = -0.5 p + 0.5 E, 

^1 

. 

1 

P + 

0 


_ e 23^ 


-0.5 


0.5 


(20) 
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Figure 3: Left and Right Laboratory Images 
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Various partial derivative estimates were computed using the expressions in (19). These values were then 
substituted in the first degree ranging polynomial to obtain the range estimates. Table 1 shows a compari- 
son between computed and actual ranges to various objects in this scene. This table also gives the number 
of points at which the ranges were computed and the standard deviation. 


It should be noted that the distribution of computed ranges corresponding to a particular object is 
not symmetric about the mean values given in Table 1. As a result, the standard deviation shown in the 
parenthesis does not truly reflect the error in estimates. A 3-D representation of the range data is shown 
in Figure 4. From this figure, it can be observed that the range estimates are more consistent in regions of 
higher contrast than in regions of low contrast. 
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Table 1: Comparison Between Actual Range (inches) and the Estimated Range (inches) 


Object 

Actual Range 

Estimated Range (Std. Dev.) 

No. of Points 

Left Pencil 

30.0 

26.99 (7.48) 

975 

Right Pencil 

26.0 

23.24 (6.33) 

2510 

Tape-on-wall 

60.0 

55.94 (8.84) 

850 

Soda can 

28.5 

27.01 (7.59) 

1738 

Bracket 

19.0 

19.71 (7.23) 

1070 



Figure 4: 3-D Representation of Computed Ranges 
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It is important to stress here that the range estimates are strongly influenced by the filter bandwidth. 
Higher bandwidth leads to partial derivative estimates that have a higher degree of fidelity. On the other 
hand, it also leads to a higher noise content at the output, often leading to inaccurate range estimates. 
Thus, as in any other signal processing application, selection of estimator parameters is largely an art if 
the noise characteristics are unknown. 

5 CONCLUSIONS 

This paper described the development and implementation of a fast technique for image-based passive 
ranging. This research was motivated by the need to automate the helicopter nap-of-the-earth flight regime. 
An optical ranging polynomial relating the scene depth with camera- vehicle parameters and various orders 
of image irradiance partial derivatives was derived. Another polynomial that predicts the truncation error 
involved in range computations was also developed. In order to enable rapid calculation of these polynomial 
coefficients, a partial derivative estimation method using causal linear filters was developed. Performance 
of the ranging algorithm was then illustrated using a laboratory image pair. 

Current research focuses on improving the accuracy of this ranging algorithm and on the inclusion of 
additional camera-vehicle motion parameters. 
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Abstract 

Nap-of-the-earth flight mode is extremely demand- 
ing on the rotorcraft pilots. This fact has motivated the 
research in automating various components of low alti- 
tude rotorcraft flight operations. Concurrent with the 
development of guidance laws, efforts are underway to 
develop systems for locating the terrain and the ob- 
stacles using inputs from passive electro-optical sensors 
such as TV cameras and infra-red imagers. A passive 
obstacle location algorithm that uses image sequences 
from cameras undergoing translational and rotational 
motion is developed. The algorithm is in a general form 
and can operate in multi-camera imaging environments. 
Performance results using an image sequence from an 
airborne camera are given. 

Introduction 

Nap-of-the-earth (NOE) flight mode is one in 
which a rotorcraft flies close to the terrain, with an alti- 
tude clearance of less than about 50 feet, w hile avoiding 
various obstacles on the way under all weather condi- 
tions. The NOE flight regime is extremely taxing on 
the rotorcraft pilots who have to carry out many other 
duties in addition to flying the aircraft. This factor has 
motivated the research in automating the NOE flight 
regime 1-3 . 

The task of locating various obstacles within the 
field-of-view has emerged as the central problem in the 
development of pilot aids for nap-of-the-earth flight. 
Use of active sensors such as Radar or laser rangers 
are not desirable due to their significant power require- 
ments and the difficulty in conducting covert missions 
with such devices on-board. This issue has driven the 
research in developing passive obstacle location schemes 
using conventional TV and infra-red imaging devices. 

Several approaches to the image-based obstacle 
location problem have emerged in recent years 4-9 . 

* Member AIAA, School of Aerospace Engineering, Georgia 
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Specifically, Reference 4 has discussed the development 
of a recursive obstacle location scheme for low altitude 
helicopter flight. In that approach, various features of 
interest in an image sequence such as regions of high 
contrast are used to determine the location of various 
objects within the field-of-view. Techniques such as 
these are termed as feature- based to distinguish them 
from techniques that do not explicitly use any object 
features for determining their location. Image-based 
obstacle location algorithms that do not explicitly em- 
ploy features are termed field-based approaches. Vari- 
ous aspects of a field-based algorithm was discussed in 
References 7-9. 

In the present work, the field-based obstacle loca- 
tion scheme discussed in References 7 - 9 is extended to 
permit the inclusion of rotational and translational mo- 
tion of the imaging devices and the rotorcraft. While 
such an extension is non-trivial, this aspect needs to be 
addressed before an operational system can be synthe- 
sized. 

Vision-Based Obstacle Location 

Passive obstacle location methods using electro- 
optical sensors have their basis in the fact that the im- 
ages obtained from sensors mounted on a moving vehi- 
cle will exhibit irradiance changes at each pixel. Rela- 
tive location of various objects, vehicle motion parame- 
ters, relative location of the imaging devices, the scene 
surface reflectance, the location of illumination sources 
all influence this irradiance change. If the surface re- 
flectance and the illumination are assumed to remain 
constant during the imaging process, then the observed 
image irradiance changes at each pixel are entirely due 
to the relative location of various objects within the 
field-of-view and the vehicle motion parameters. Addi- 
tionally, if all the objects within the field of view are 
assumed fixed with respect to an inertial frame, then 
the irradiance change at each pixel location can be di- 
rectly related to the location of these objects. In this 
case, if the vehicle motion parameters and the imaging 
device constants are known, it is possible to determine 
the location of various points within the field of view. 

Image-based obstacle location algorithms operate 
on the basis of the correspondence hypothesis . The 
central idea here is to establish the correspondence be- 
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tween various objects in a pair of images and to measure 
the displacement of these objects on the image plane. 
The measured displacement or the disparity can then 
be employed together with the perspective projection 
geometry for computing the location of various objects 
within the field-of-view. 

In order further to crystallize these ideas, consider 
the various coordinate systems illustrated in Figure 1. 
The first of these is the image plane coordinate system 


Ob>*ct 



with the major axis of the image plane being desig- 
nated as the X p axis and the minor axis being labelled 
the Y p axis. The origin of this image plane coordinate 
system is located at the center of the aperture. In this 
coordinate system, an image may be defined as an ir- 
radiance distribution E(x pf y p ) ) with E being the ir- 
radiance specified on a gray scale at a point x p ,y p on 
the image plane. Next consider the camera coordinate 
system X Cy Y C} Z c with the X c axis passing through 
the origin of the image-plane coordinate system, and 
the axes Y c , Z c being parallel to the image plane co- 
ordinate system. The origin of the camera coordinate 
system is at the lens center, located one focal length / 
behind the image plane X p , Y p . Such an arrangement 
of the camera coordinate system is routinely employed 
in image processing work to avoid having to deal with 
inverted images. 

The origin of the camera coordinate system is as- 
sumed to be located at a point l \ , / 2 , I 3 with respect to 
the body axis, and oriented by three angles e 1} c 2 , £3 
about the Z Cy Y c and X c axes respectively. The vehicle 
body axis system is defined using the standard flight 
dynamics convention, viz., the X * axis pointing along 
the nose of the rotorcraft, Y axis along the starboard 


direction, and the Zb axis completing the right handed 
triad. The body axis system may be related to an in- 
ertial frame X , Y, Z through yaw, pitch, roll Euler 
angles 9 , <j> and the three translational components 

x viVvyZv The definition of the inertial frame follows 
the standard flight dynamics convention with the X- 
axis pointing towards north, the Y-axis pointing to- 
wards east and the Z-axis pointed in the direction of 
local gravity vector. 

A point x py y p corresponding to the object in the 
image plane is related to its location with respect to 
the inertial frame x,y, z through various translations 
and rotations of the vehicle body frame and the camera 
frame, together with the perspective projection. 

Let E\(x p ,y p ) and E 2 (x py y p ) be two different 
views of a sample scene. Assuming that these two im- 
ages contain all the objects of interest and that to a 
large degree, the perceived irradiance of the scene has 
remained invariant during the imaging process, the cor- 
respondence hypothesis can be expressed as 

E\(x p ,y p ) = E 2 (x p + Ax p ,y p + Ay p ) (1) 


E 2 (x p ,y p ) = Ei{x p - Ax p ,y p - Ay p ) (2) 

Here, Ax p ,Ay p are the disparities between the corre- 
sponding points in the two images defined at every point 
on the image plane. 

It is important to stress that the correspondence 
hypothesis expressed in (1) and (2) should be inter- 
preted in a limited sense because it does not account 
for various occlusions that can arise during the imaging 
process. These occlusions arise due to the the viewing 
angle and the camera aperture dimensions. Note that 
equations such as (1) and (2) can be written for any 
number of images, provided that the objects of interest 
appear in all the images. Various methods for satisfy- 
ing the correspondence hypothesis have been discussed 
in the literature 4-10 . In the present research, the cor- 
respondence hypothesis will be approximated by first 
expanding the expressions (1) and (2) in a Taylor se- 
ries and then truncating them based on the acceptable 
computational complexity. Such an expansion treats 
disparities Ax py A y p as perturbations about the spec- 
ified pixel locations. Clearly, the error in such an ap- 
proximation depends on the number of terms included 
in the Taylor series. 

In all that follows, a bold face lower case letter 
will denote a vector, while a bold face upper case letter 
will denote a matrix. A superscript T will be used to 
denote the vector-matrix transpose operation. Define 
the disparity vector d, the irradiance gradient vector g, 
and the irradiance Hessian matrix H as 


A x p 


' dE/dxp ' 

. A y P 

i 8 — 

. dE/dy p 
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d 2 E/dx p 2 d 2 Efdx p dy p 

d 2 E/dy v dx p d 2 E/dy p 2 


(3) 


Now, if are the position vector components of 

the point x c , y c , z c in the body axis system, 


Equations (1) and (2) can now be expanded using 


a two-dimensional Taylor series to yield 
E\ — E 2 — g2<l + — d^H2d + 

Similarly, 


X c 


Xb i 1 

Vc 

= Tj 

yk - h 

m 


Zi-l 3 


where 


( 10 ) 


Ei-Ei = -g 1 d+ld T H 1 d- (5) 

Expressions (4) and (5) relate the irradiance spatial 
partial derivatives to the irradiance change and the im- 
age disparities at every pixel. It is awkward to express 
the terms beyond second-order in (4) and (5) using the 
vector-matrix notation because they contain Tensors. 

In the special case where the disparities Ax pi Ay p 
arise from the vehicle motion resulting in apparent im- 
age velocities, the equations (4) and (5) with the first 
term on the right hand side is the optical flow constraint 
equation 10 . In this case, the apparent image velocity 
components are called the optical flow velocities. This 
characterization is avoided here because the disparities 
Ax p , Ay p may arise from other sources also, stereo 
camera arrangement being an example. By adopting 
the point-of-view proposed in the foregoing, equations 
(4) and (5) may be used for temporal image sequences 
from a single camera motion sequence as well as simul- 
taneous images obtained from several cameras. 

Adding and subtracting the expressions (4) and (5) 
results in 


Ti - 

CC2CCI CC^SC 1 —SC 2 

— CC 3 SC 1 SC 3 Se 2 Sei -f Cl 3 Cl\ SC3C62 
CC3SC2CCI + sessti ce3se 2 sei - se 3 cei ce 3 cc2 

(ii) 

The variables c and s in the 3x3 matrix in (11) are the 
sine and cosine functions. In a more concise notation, 

x c = T 1 (x6-1) (12) 

Next, if the location of the same object with re- 
spect to an inertial frame is x, y, z, and the vehicle is 
located at a point x Vi y v , z v with respect to this iner- 
tial frame, the position components of the object point 
in the body frame are given by 


" x b " 


’ X - x v 

Vb 

= t 2 

y-yv 

z *> 


z — z* 


where 







c0c\!) 

cOsxj) 

-s0 

r i, i lT 

■ 

d+ .... 

... (6) 

t 2 = 

$<f>s9cip — c<f>srp 

s(f>s9s 4- c4>cip 

s<f>c6 

gl+g2 d+jd 

h 2 -h, 


c4>$0cxl> + S<f>S\j) 

c<f>s0sip — S<f)Ctp 

c<j>c6 


°= g2-gijd+ld 7 


Hi + H 2 


d + 


(7) 


(14) 

Here, 0, <t> are the yaw, pitch, and roll attitudes of the 
rotorcraft. Equation (13) can be written in a compact 
form as 


Next, as in Reference 8, equation (6) will be used 
for obstacle location, while the expression (7) will be 
used for computing the truncation error in the Taylor 
series approximation. The components of the disparity 
vector Ax p , Ay p can now be related to the vehicle mo- 
tion parameters and the location of various objects in 
the field-of-view. 

Let x c , y c , z c be the location of a point on an object 
with respect to the camera coordinate system. Since 
origin of the image plane is located at the point [/ 0 0] T 
with respect to the camera frame, this object point 
would appear at a point x p ,y p on the image plane. If 
the camera focal length is /, then perspective projec- 
tion rules require that 


X6 = T 2 (x-x t ,) (15) 

Since the objective is to eliminate the disparity vec- 
tor d in favor of the object position vector components, 
equations (6) and (7) examine the changes in x p , y p 
in response to the changes in the vehicle location by 
Ax v , Ay v , A z v . For the sake of clarity, this can be 
carried out in two steps. First consider the changes in 
x p , y p in response to changes in x c , y C} z c . This can 
be accomplished by evaluating the expressions (8) and 
(9) with the camera referenced object location being de- 
fined as x c + Ax c , y c +Ay c , z c +Az c . Using elementary 
algebraic operations it can be shown that: 


Xp _ yc 
/ x c 

(8) 

Ax p = f 

yp _ £ £ 
/ x c 

(9) 

II 

<1 


Ay c 

Xc -h Ax c 
Az c 

x c + Ax c 


Xn 


Aj c 

x c + Ax c 


Ax c 

Vp x c + Ax c 


(16) 

(17) 
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These changes are assumed to occur solely due to 
the motion of the rotorcraft. Next examine the changes 
in z c , y Cl z c due to the changes in vehicle position com- 
ponents x v , y v , z v . Using the coordinate transforma- 
tions discussed (12) and (15) one has that: 

Ax c = TiAxj (18) 

Ax* = AT 2 (x - x„) - T 2 Ax, - AT 2 Ax v (19) 

The matrix AT 2 is the change in the transforma- 
tion matrix relating the body frame and the inertial 
frame caused by the changes in the vehicle attitude. 
Note that the equation (18) assumes that the camera 
is fixed relative to the vehicle . Equation (19) assumes 
that the obstacles are fixed with respect to the inertial 
frame. Next substituting (19) in (18) results in: 


Note that the Ax p , A y p are bilinear with respect 
to the inverse of the object position component x c mea- 
sured with respect to the camera coordinate system. 
Next, the expressions (25) and (26) may be substituted 
in the Taylor series approximation of the correspon- 
dence equation (6) to yield an equation that relates the 
vehicle-camera motion parameters with the image ir- 
radiance partial derivatives. In order to simplify the 
development, only a linear approximation will be dis- 
cussed in this paper. Approach to include second-order 
terms is direct. From (6), the first-order Taylor series 
approximation of the correspondence hypothesis is: 


1 

dE 1 

dE 2 

A 1 
A /yi 1 _ 

'dE 1 

oe 2 

2 

dx p 


Ax p + 2 

,dy p 

dVp. 


Substituting for the disparity components Ax py 
A y p in expression (27) from (25), (26) yields: 


Ax c = T 1 [AT 2 (x - x„) - (T 2 + AT 2 )Ax,] (20) 
Using (12) and (15), 

X - x„ = T 2 T [T 1 t x c + 1] (21) 

Substituting equation (21) in (20) yields: 


Ax c = T 1 [AT 2 T 2 t (T 1 t x c + 1) - (T 2 + AT 2 )Ax„] 

( 22 ) 

Finally, dividing expression (22) by x c on both 
sides: 


Ax c /x c ' 


1 

Ay c /x c 

= T 1 AT 2 T 2 t Ti T 

Xp/f 

Az c /x c 


. yp/f . 


+T,AT 2 T 2 t 


h/x c 


A 2 ?^ 1 x c 

h/x c 

-T,(T 2 +AT 2 ) 

A y v /x c 

H 


Az v j x c 


(23) 

For the sake of brevity, the equation (23) can be 
expressed in the form 


_ *3^1 ~ ai(/*4 ~ x pk 3) ~ 6i(/fc 5 ~ Vp k 3 ) 

ai(fki - x p k 0 ) 4 bi(fk 7 - y p k 0 ) - ci(l + k 0 ) 

(28) 

where 


1 

'dE 1 dE 2 

2 

dx p + dx p 

1 

dEi 8E 2 

2 

. 9y P + dy p . 


(29) 

(30) 


c\ = Ei — E\ (31) 

Various quantities on the right hand side of equa- 
tion (28) can be obtained from on-board instruments 
and the given image sequences. Hence, equation (28) 
can be solved to obtain x c . Once x c is computed at ev- 
ery pixel location, the two remaining components of the 
object position vector at every pixel can be determined 
using the relations: 


y c = x c x p /f , z c = x c y p /f (32) 

The position of the object with respect to the in- 
ertial frame is then given by: 


' Ax c /x c ‘ 


ko 4 k 3 /x c 

Ay c /x c 


*1 4 k 4 /x c 

. Az e /x c _ 


^2 4 £ 5 / x c 


(24 


with parameters fco, £ 1 , &2> ^ 3 > ^5 being computed us 

ing the operations indicated in equation (23). The row: 
of equation (24) can now be substituted in the expres 
sions (16) and (17) to yield 


A _ fkj Xpkp 4 (/&4 X p k 3 )/x c 

P ~ 1 4 ko 4 k 3 /x c (25) 


* _ fkz y P k 0 + (fk 5 y P k 3 )/x c 

- r+to + 4 3/x c (26) 


X = x„ 4 T 2 T (Ti T x c -I- 1) (33) 

The instantaneous vehicle position vector Xt,, the 
instantaneous transformation matrix from the vehicle 
body frame to the inertial frame T 2 , its change AT 2 , 
the constant transformation matrix from the vehicle 
body frame to the camera frame Ti, the camera fo- 
cal length / and its orientation angles are known from 
various on-board instruments. 

The image sequence related quantities in the ex- 
pression (28) are the image spatial partial derivatives 
dEi/dx p , dE\/dy p) dE 2 /dx p) dE^/dyp and the irradi- 
ance difference £ 2 - E\ at every pixel. A method for 
estimating these quantities will be presented in the next 
section. 


47 



Partial Derivative Estimation 

Several numerical techniques are available in 
the literature for the computation of the partial 
derivatives 10 . These range from simple central differ- 
ence schemes to very sophisticated algorithms that re- 
ject noise and unwanted frequency components in the 
input signal. It is well known in the signal processing 
literature that derivative estimation is a noise amplify- 
ing process. Given an irradiance distribution E(x py y p ), 
the objective of the partial derivative estimators are 
to provide sufficiently noise free estimates of dE/dx p , 
dE/dy Pi d 2 Ejdx 2 y d 2 E/dy p 2 , d 2 E/dx p dy p , .... to- 
gether with a consistent image irradiance distribution. 
Additionally, it is desirable to carryout this estimation 
process as fast as possible to enable real-time or near 
real-time implementation of the object location algo- 
rithm. 

Frequency domain estimation methods are popular 
in the image processing literature 10 . While these ap- 
proaches enable a direct formulation of the non-causal 
estimation problem, they are unattractive for use with 
electro-optical sensors that form images using a sequen- 
tial scanning process. In the following, a partial deriva- 
tive estimator will be formulated directly in the spa- 
tial domain. Unlike the frequency domain estimators, 
the present formulation will employ causal estimators 
to permit sequential computations. Although this ap- 
proximation compromises accuracy to a certain extent, 
numerical experiments have shown that the error mag- 
nitude can be controlled by an appropriate selection of 
the estimator parameters. 

Following the frequency domain image processing 
methods 10 , the image irradiance distribution E(x p ,y p ) 
will be assumed to be factored into a product of two 
functions of one variable each, denoted by the sym- 
bols ei(xp) and e 2 (t/p )- This factorization allows the 
evaluation of double integrals as products of two, sin- 
gle integrals. Note that such a factorization is routinely 
employed to produce solutions of linear partial differen- 
tial equations 11 . In the following, estimated quantities 
will be denoted by a ? Thus, 


With this factorization, the estimation problem 
can be cast as two, coupled, one-dimensional estima- 
tion problems. Although the estimation problem can 
be formulated either as a serial process or as a parallel 
process 9 , a serial implementation will be discussed in 
this paper. 

Define p as the n dimensional estimator state vec- 
tor. The estimator dimension n may be selected based 
on the order of partial derivatives to be estimated and 
the desired smoothness requirement. For instance, if 
it is desired to estimate the first-order and the second- 
order partial derivatives with the requirement that the 
second-order partial derivatives be smooth, n has to 
be at least 3. Next consider a linear dynamic process 
evolving along the x p direction with the measured im- 
age irradiance E as the input. 

— p = A p + B E (36) 

This dynamic process is assumed to evolve along 
the image scan lines. The matrices A and B can be cho- 
sen to obtain the desired transmission gain and band- 
width for the x p spatial frequencies. The output equa- 
tion for this linear dynamic process can be defined as: 

r=C P + D £ , r = «,fo ^0 f (37) 

The matrices C, D should be chosen based on the lin- 
ear dynamic system (36) dimension and on the desired 
orders of the partial derivatives. 

Similarly, define a linear dynamic process evolving 
along the y p direction, with the first element of the 
vector r as the input, i. e., 

d 

— q = M q -f N e 2 ei (38) 

The matrices M and N can be chosen to obtain the 
desired gain and bandwidth for the y p spatial frequency 
components. The output equation of this dynamic pro- 
cess can be chosen as 


E= ei(xp) e 2 (y p ), E = ei(x p ) e 2 (y p ) (34) 

E is the estimated image irradiance distribution. The 
partial derivatives of E can be written in terms of the 
total derivatives of its factors e\, e? as: 


dE__ . de L dE_ _ „ de? 

dx p 2 dx p dy p 61 dy p ’ 

d*_E _ . d 2 E _ _ 

dxp* _C2 dxp 2’ dt -£l dyp 7 ' 

d\E _ de i de 2 
dx p dy p dxp dy p ' 


s = Kq + Le ie2 , s = e l [e 2 p-pl ,f (39) 

dy P dy p 2 

with the matrices K, L being chosen in the same man- 
ner as the output matrices C, D in equation (37). The 
vectors r and s can now be used to compute the esti- 
mated image irradiance and its various partial deriva- 
tives. Elements of these vectors will be denoted by lower 
case letters in the following. It can be verified that 



S\t2 dE ^ 
r i ' dyp 52 ’ 


d 2 E _ $ir 3 d 2 E __ 

dxp 2 n ’ dtf ~ S3 
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d 2 E r 2 s 2 
dx p d y p r\ 


(40) 


In equation (40), it is important to ensure that the 
first element of the r vector r\ does not become zero at 
any point on the image. This can be assured by adding 
a constant bias to the input image E and subtracting 
this bias at the output. Note that such an operation 
does not alter the partial derivative estimates. 

Since the partial derivative estimation process was 
formulated keeping the image forming process in view, 
these quantities are available as the image is being 
formed. As a result, by the time the image scan- 
ning process is complete, various parameters required in 
equation (28) can be computed. Thus, after the avail- 
ability of the first frame, the object location computa- 
tions can be completed at the camera frame rate. 

A block diagram of the object location algorithm 
is given in Figure 2. This block diagram assumes that 



Figure 2: The Object Location Algorithm 

an image pair is available at the beginning of the com- 
putations. The estimators along the x p and y r direc- 
tions provide the noise attenuated partial derivatives to- 
gether with consistent image irradiance estimates. The 
coefficients of equation (28) are then computed by com- 
bining this information with the vehicle-camera motion 
data. The position vector component x c corresponding 
to every point on the image plane is then computed by 
applying the expression (28) at every pixel. The ve- 
hicle and the camera rotational and translational data 
required for this calculation is assumed to be available 
from on-board instruments. The other components of 
the position vector y c , 2 C are then determined using the 
expression (32) at every pixel. 


The quantities x Cl y Cl z c are computed only at 
regions where the inter-frame irradiance changes and 
the spatial partial derivatives are sufficiently large. The 
remaining points are assigned a pre-defined maximum 
scene depth. Once the position vectors in the camera 
frame become available, the expression (33) can be used 
to transform the position components into the inertial 
frame. 


Results and Discussions 

The performance of the obstacle location algorithm 
was tested using two different image sequences. The 
first of these was a laboratory scene used in References 
8 and 9. The second test image sequence was obtained 
from a camera mounted on the nose of a helicopter. 
Results obtained using the second image sequence will 
be discussed in this paper. Two images from the air- 
borne camera are shown in Figures 3 and 4. These 



Figure 3: The First Sample Image 

images consist of a 512 x 512 pixel arrays, with 8 bit 
gray-scale digitization. The camera was operating at a 
rate of 30 frames/second. During the imaging process, 
the rotorcraft was flying at an altitude of about 12 feet 
above the runway at a speed of 32.6 feet/second. The 
images in Figure 3 and 4 are temporally separated by 
0.17 seconds. During this time, the rotorcraft under- 
went a translational motion of [5.41, —0.44, 0.02] T feet 
and experienced an attitude motion of 0.059 degrees 
about the pitch axis, -0.08 degrees about the yaw axis, 
and 0.12 degrees about the roll axis. 

Besides other things, the Figures 3 and 4 show a 
runway, and five vehicles parked on the two sides. The 
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Figure 4: The Second Sample Image 


Figure 5: Points at which Object Location data is Com- 
puted 


instrumented nose boom of the helicopter is clearly vis- 
ible at the top right hand part of the image. Time 
sync, parameters identifying this frame are displayed 
in a sub-picture on the top left-hand-side of the image. 

In order to compute the spatial partial derivatives 
together with a consistent irradiance distribution, these 
images were processed though irradiance estimators of 
the form: 


-T- = —5.0 p + 5.0 E } 

dXp 



j — — —5.0 p + 5.0 Cje2, 

<*Vp 



Output of these estimators were used in equation 
(40) to form the various partial derivatives. These val- 
ues were then used together with the vehicle motion 
data obtained from the on-board inertial navigation sys- 
tem to compute the distance to various obstacles. 

The object location data was computed at 480 
points on the image plane. In Figure 5, the points at 
which the location data was computed are shown as 
white dots. This represents about 0.2 % of the possible 
262,144 points. The object location data can be reli- 
ably recovered only at those points where a significant 
image intensity change is detected between frames. 


Object 

Actual 

Range 

Estimated 

Range 

Number 
of Points 

Vehicle 1 

356.4 

335.6 

113 

Vehicle 2 

607.0 

395.9 

24 

Vehicle 3 

232.8 

202.1 

160 


Table 1. Comparison between Actual Range (feet) and 
Estimated Range (feet) for Some of the Objects in the 
Scene 

Table 1 summarizes the comparison between com- 
puted and actual locations of some of the objects in 
this scene. This table also gives the number of points 
at which the computations were carried out for each 
object. Eventhough the image was not of high qual- 
ity, note that the algorithm determines the location of 
two vehicles with acceptable accuracy. There is a large 
error in the location of the vehicle number 2. At this 
stage it appears that this inaccuracy may be caused 
by the low resolution of the image. Additional tests 
are currently being carried out to determine the image 
resolution requirement for obtaining a desired position 
determination accuracy. 


Conclusions 

A passive image-based obstacle location algorithm 
for use in the guidance of rotorcraft during nap-of-the- 
earth flight was described. The algorithm uses the 
vehicle-camera translation and rotational motion data 
together with the image sequences to compute the lo- 
cation of various objects within the field of view. 
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The method uses a Taylor series approximation of 
the correspondence hypothesis together with the per- 
spective projection geometry. The present research gen- 
eralizes previous work on an object location algorithm 
to include both the translational and rotational motion 
of the imaging devices and the vehicle. The resulting 
algorithm is more general and can operate in multi- 
camera imaging environments. The performance results 
using an image sequence from an airborne camera are 
given. 
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5. Derivative-Free Vision-Based Ranging 

In all the computational schemes discussed in this report, the parameters enabling ranging were 
the spatial partial derivatives of the image. It is well known in the signal processing discipline that 
derivative estimation is a noise amplifying process. Consequently, if a derivative free method can be 
synthesized, then it would be less sensitive to picture noise. 

This idea forms the basis for the ranging scheme discussed in the following paper. In this scheme, 
instead of satisfying the correpondence hypothesis using partial derivative estimates as in previous 
chapters, the integral error between a predicted value of image irradiance and those from the actual 
measurements are minimized. The solution is required to satisfy a differential constraint derived by 
approximating the correspondence hypothesis using Pade expansion. 

Minimization of the integral irradiance error subject to the differential constraint results in a 
calculus of variations problem. The necessary conditions for a minimum yield a set of linear two-point 
boundary-value problems, which can be solved using the backward sweep method. In this process, 
the sum and the difference between two images are used to form a Riccati equation and a linear 
differential equation. These equations are integrated from the right edge of the image to the left edge 
to compute the feedback gains for the range estimator. Multiplying these feedback gains with the 
difference between the actual and the predicted image irradiances then produce the range estimates. 

A paper outlining this work has been communicated for presentation at the 1992 AIAA Guidance , 
Navigation, and Control Conference , and is included in the following. 
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Vision-Based Stereo Ranging 
as an Optimal Control Problem 
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Abstract 

The recent interest in the use of machine vi- 
sion for flight vehicle guidance is motivated by 
the need to automate the nap-of-the-earth flight 
regime of helicopters. Vision- based stereo rang- 
ing problem is cast as an optimal control problem 
in this paper. A quadratic performance index con- 
sisting of the integral of the error between observed 
image irradiances and those predicted by a Pade’ 
approximation of the correpondence hypothesis is 
then used to define an optimization problem. The 
necessary conditions for optimality yield a set of 
linear two-point boundary-value problems. These 
two-point boundary- value problems are solved in 
feedback form using a version of the backward 
sweep method. Application of the ranging algo- 
rithm is illustrated using a laboratory image pair. 

Introduction 

Passive ranging has emerged as a central issue 
in nap-of-the-earth helicopter flight guidance 1 " 3 
and autonomous planetary rover mission 4 " 6 as ev- 
idenced by the recent appearence of several papers 
in this area. Research on both the sensors and 
guidance algorithms are currently underway. The 
status of research activity on the nap-of-the-earth 
guidance problem is summarized in Reference 3. 
The research reported in this paper is a compo- 
nent of the helicopter guidance problem. 

Image-based passive ranging is an active re- 
search area in the robotics discipline' . Applica- 

* Member AIAA, School of Aerospace Engineering, Geor- 
gia Institute of technology, Atlanta. Mailing Address : FSN 
Branch, M.S. 210-9, NASA Ames Research Center 

* Research Scientist, Member AIAA, FSN Branch, MS 
210-9 

♦Member AIAA, Member of the Professional Staff, Ster- 
ling Software, Inc. 


tion of this technology for flight vehicle guidance 
is a more recent development. Driven by the need 
to determine range to various obstacles within the 
field-of-view, Reference 8 outlines the development 
of an image-based recursive range determination 
scheme for nap-of-the-earth helicopter flight guid- 
ance. In that approach, various features such as 
edges or regions of high contrast in an image se- 
quence were used to determine the range to var- 
ious objects within the field-of-view. Techniques 
such as this are termed as feature-based ranging to 
distinguish them from techniques that do not ex- 
plicitly use any features for range determination. 
Field-based ranging algorithms, on the other hand 
do not explicitly employ any features for range 
computations. Various aspects of the development 
of a field-based ranging algorithm was described 
in References 9 - 12. The algorithm developed in 
these papers relied on the image irradiance par- 
tial derivatives to compute the range. The anal- 
ysis considered both motion and stereo image se- 
quences. 

Present paper reports on the development of 
a field based stereo ranging scheme that does not 
require spatial partial derivative estimates. This 
research develops a differential constraint based on 
the approximation of the correspondence hypoth- 
esis using Pade’ expansion in a registered stereo 
pair. Under certain conditions, this expansion can 
be shown to be identical to the Taylor series expan- 
sion employed in References 9 - 12. Using perspec- 
tive projection relationships, the image disparity 
is next eliminated in favor of the stereo baseline, 
camera focal length and range. The state variable 
in the resulting differential constraint is the sum 
of the image irradiances and the control variable 
is the range to various points within the field-of- 
view. Integral of a quadratic form in the error 
between the measured irradiances and those pre- 
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dieted by the Pade’ approximation of the corre- 
pondence hypothesis is then defined as the perfor- 
mance index. Necessary conditions for minimizing 
the performance index subject to the differential 
constraint are obtained by the application of opti- 
mal control theory 13 . This process produces a set 
of linear two-point boundary- value problems which 
can be solved for range in state feedback form us- 
ing a version of the backward sweep algorithm 13 . 

The ensuing sections will discuss the develop- 
ment of this algorithm in further detail. The per- 
formance of the optimal stereo algorithm will be 
demonstrated using a laboratory image pair. 

Stereo Ranging 

Images obtained from a pair of spatially sep- 
arated cameras will exhibit irradiance differences 
at each pixel. These differences depend on t ho rel- 
ative location of various objects, relative location 
of the cameras, the scene surface reflectance and 
the location of illumination. If the scene surface 
reflectance is relatively invariant between the two 
images, the reason for observed image irradiance 
changes at each pixel may be attributed to the rel- 
ative location of various objects within the field-of- 
view. This fact forms the basis for correspondence 
hypothesis. 

According to the correspondence hypothesis, 
the relative displacement of the corresponding ob- 
jects in an image pair or the disparity can be em- 
ployed together with the perspective projection 
geometry to compute the range to these objects. 
Following the existing literature, the frontal im 
age plane representation will be used for image 
description. The ensuing analysis will deal with 
monochromatic stereo pairs obtained using identi- 
cal cameras. 

Consider a registered stereo pair together with 
the coordinate systems as shown in Figure 1. In 
Figure 1, the major axes of the two image pianos 
are designated as the X p axis and the minor axes 
are labelled the Y p axis. The origin of this image 
plane coordinate system is located at the center of 
the aperture. The range ^ to various objects within 
the field-of-view will be measured with respect to 
a coordinate system fixed to the origin of the first 
camera lens, located one focal length / behind the 
image plane with the X and Y axis being parallel 
to X p and Y p axes. The second camera is offset 



Figure 1: The coordinate system. 


from the orgin of this sytem by a distance Ax along 
the X axis. The distance Ax is termed as the 
stereo baseline. 

In the image plane coordinate system, an im- 
age may be defined as an irradiance distribution 
E(x p , Up) in the X p , Y p plane, with E being the ir- 
radiance specified on a gray scale at a point (x p , y p ) 
on the image plane. It is assumed here that except 
at finite number of points, the image irradiance 
distribution is continuous along the horizontal di- 
rection x p . This is a reasonable assumption be- 
cause imaging sensors perform spatial integration. 
Since the analysis assumes registered stereo pairs, 
any number of discontinuities are permitted along 
the y p direction. For the purposes of the present 
analysis, nonlinear effects contributed by gray scale 
quantization will be ignored. 

Let E\{x p , y p ) and E 2 {x p , y p ) be a given stereo 
pair obtained from cameras with identical param- 
eters. Assume that the two images contain all the 
objects of interest and that to a large degree, the 
irradiance of the scene is independent of the cam- 
era position. Further assume that the two images 
are perfectly registered 7 . In this case, the corre- 
spondence hypothesis may be expressed as 

Ei(x p ) - E 2 (x p + A x p ), y v - constant (1) 

This equation states that an object at the point 
{x p ,y p ) in the first image will appear at the point 
(x p -f A x p ,y p ) in the second image. Here, Ax p is 
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the disparity between the corresponding objects in 
the two images. Since different points in the image 
plane may correspond to different objects within 
the field-of-view, the disparity Ax p may vary as a 
function of x v . The depth information is encoded 
in these disparities. Note that the disparities in 
the y p direction are zero because the stereo pair is 
assumed to be registered. Further, note that the 
limited camera field-of-view and the various occlu- 
sions that can arise in real scenes may result in the 
correspondence hypothesis not being valid at every 
point on the image pair. 

Pade’ Approximation 

The correpondence hypothesis may be satisfied 
in various ways. Indeed, the method employed 
for satisfying the correspondance hypothesis is the 
main feature that distinguishes one stereo ranging 
algorithm from another. The feature-based family 
of methods 8 attempt to decrease the complexity 
involved in satisfying the correspondence hypoth- 
esis by first extracting the features of interest in 
an image such as edges or regions of high contrast 
and then establishing the correspondence between 
them. This step reduces the correspondence task 
by a significant amount. Pixel by pixel correlation 
methods have also been proposed for satisfying the 
correspondence hypothesis 14 . 

In the present paper, an alternate approach 
for satisfying the correspondence hypothesis will 
be advanced. Using Laplace transforms 15 together 
with the shift theorem, the correspondence equa- 
tion (1) over small patches can be written as: 

Ei(s) = e sAx »E 2 (s) (2) 

Here s is the Laplacian operator. According to 
equations (1), and (2), the image irradiance distri- 
bution in the second image is simply the image ir- 
radiance in the first image shifted by the unknown 
disparity at a pixel. Expression (2) can next be 
expanded using Pade’ approximation 16 . The ex- 
pression (2) is first rewritten as 

sAip Aip 

e-^ E E i (s) = e s ~^E 2 (s) (3) 

Expanding both sides, 

Ei(s) - ^sEi{s)Ax p + 

= E 2 {s)+ l -sE 2 (s) Ax p + (4) 


Retaining the two leading terms on both sides and 
inverse transforming yields the differential equa- 
tion 


2 (E 1 -E 2 )-±- = -}^(E 1 + E 2 ) (5) 

A similar expression was derived in Reference 
10 using forward and backward expansion of the 
Taylor series. In order to assure continuous depen- 
dence of the left hand side of the expression (5) on 
Ax p , it will be assumed that E\ - E 2 7^ 0 over 
any finite interval of x p . It is sometimes useful 
to retain higher-order terms in the Pade’ expan- 
sion. However, the resulting differential equation 
for correspondence hypothesis will be nonlinear in 
1/A.Tp, making the solution process more complex. 

Next, the geometry of image forming process 
will be used to eliminate the unknown disparity 
A x p from expression (5) in favor of camera param- 
eters and range. From the perspective projection 
geometry, an object located at a position (,r, 2) 
will appear at a point 


Xp 

/ 


x 


(6) 


in the first image. Since the second image in the 
registered stereo pair is located at (x + A.t), the 
same object will appear at a point 


x p + A x p x + Aa; 

7 


(7) 


in the second image. Ax p is the disparity at a 
point on the image plane corresponding to an ob- 
ject. An expression relating the disparity, stereo 
baseline and the range to an object can be obtained 
by subtracting expression (6) from expression (7), 
Viz. 


A* p = /— (8) 

Note that the stereo baseline Ax and the cam- 
era focal length / are known. Expression (8) re- 
lates the observed disparity Ax p at a pixel location 
to the range to an object z. Using expression (8) 
in the approximation for the correspondence hy- 
pothesis (5) yields 


-E(Ey + E 2 ) 

(IX p 


2 [E x 


/Az 

Next, define two new variables 


(9) 
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G = E x + E 2 , /I = (10) 

JLiX 

So that the expression ( 9 ) becomes 

c; (it) 

From this point on, a dot over variables will be 
used to denote differentiation with respect to the 
independent variable x p . Expression (11) relates 
the sum and differences of the observed image ir- 
radiances with the range to various objects within 
the field-of-view. If the image irradiances, focal 
length and the stereo baseline are given, then the 
coefficient A(x p ) can be computed at every point 
on the image plane. In this case, expression (11) 
is a varying coefficient first-order linear differential 
equation with z as the input and G as the state. 

In the next section, the differential constraint 
expressed by equation (11) will be used together 
with the measured image irradiance to determine 
the range as a function of x p . 

Ranging as a Two-Point Boundary- Value 
Problem 

The data available in a registered stereo pair are 
the image irradiances E\{x p , y p ) and E 2 (x p , y p ) 
along every y p within the camera frame. Accord- 
ing to the differential constraint approximating the 
correspondence hypothesis given (11), at every y p 
coordinate, the derivative of the sum of the ir- 
radiances when divided by the varying coefficient 
A{x p ) will yield the range 2 as a function of the 
image plane coordinate x p . References 9-12 ap- 
proched the ranging problem in precisely this man- 
ner. Specifically, in Reference 12, the image irra- 
diance derivatives were first computed using linear 
estimators. An approximation to the correspon- 
dence hypothesis such as equation (11) was then 
used to compute range. In the following, an alter- 
native to this two step procedure will be developed. 

Since the differential constraint (11) is only an 
approximation and because of the fact that the ob- 
served irradiances are corrupted by noise, one may 
require that this differential constraint be satisfied 
in a minimum integral square error fashion. To 
this end, consider the problem 


min 




I xp=x p (r) 



subject to the differential constraint 

G = A z ( 13 ) 

G is the sum of the measured irradiances, and the 
coefficient A(a:p) is computed from the measured 
irradiances and camera parameters using the ex- 
pression (10). The upper and lower integration 
limits x p (l), x p { r) correspond to the left and right 
image boundaries. Although by no means essen- 
tial, in order to simplify the ensuing development, 
it will be assumed that the sum of the measured 
irradiances G and G computed from the differen- 
tial constraint ( 13 ) are equal at the left boundary 
of the image. This implies that the left boundary 
condition G(x p (l)) = G(x p (l)) is specified, while 
the right boundary condition is free. 

The first term weighted by the parameter 7 
in (12) helps to ensure that the measured sum of 
irradiances G and the predicted sum of irradiances 
G are close at the right edge of the image. The 
weighting factors a and (3 in the performance index 
(12) can be used to establish the trade-off between 
irradiance error and the range. Bryson’s rule 13 can 
be used to select the weighting factors as : 


a S 


1 

( 14 ) 

[ X p( T ) 

- x p {1))Gm 

0 ^ 


1 

( 15 ) 

[x P (r) 

- X p (1)]zm 

7 


1 

( 16 ) 

Gm 

\x p =x p (r) 


where Gm is the maximum acceptable value of 
(G — G ) 2 and z\f is the maximum acceptable value 
of z 2 . 

Attention is drawn to the fact that the second 
term in integrand of (12) is essential to ensure the 
existence of a non- trivial solution. This term has 
the effect of limiting the magnitude of the range 
with respect to irradiance changes. The weighting 
factors q, (3 are permitted to be functions of x p . 
It is assumed here that the computation proceeds 
from the left boundary of the image to the right 
boundary. 

The expressions (12) and ( 13 ) define an opti- 
mal control problem, with G as the state variable 
and 2 as the control variable. Necessary condi- 
tions for a minimum can be obtained by formally 
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applying Pontriagin’s minimum principle 13 . The 
first step here is the formulation of the variational 
Hamiltonian 13 as : 

H = ~[G-C] 2 + ^~ 2 + XA z (17) 

A is the unknown costate. The necessary condi- 
tions for optimum are 13 : 

A = -a[G -G] (18) 


XA 

T 


(19) 


Since the sum of the irradiances G is required 
to satisfiy the measured value G exactly at the left 
boundary of the image plane, the costate A at this 
boundary is unspecified. At the right boundary, 
the term minimizing the square of the deviations 
on the right boundary dictates the value of the 
costate, i. e. 


A(z p (r)) = 7 [G - G] \ Xp=Xp{r) (20) 

Substituting for £ in ( 13) using expression ( 19) 
yields the two-point boundary- value problem at 
each specified value of y p as : 

A 2 

G = - — A, G{Xp(l)) Specified (21) 
A = -a[G - G\, 

A (x p (l)) unknown , A(x p (r)) = 7(6' - G]\ x = x ^ r ) 

( 22 ) 

The linear two-point, boundary value problem 
(21), (22) can be solved in feedback form for G and 
A by using a version of the backward sweep 13 , as 
will be demonstrated in the following. Substituting 
these values in expression (19) then yields range at 
each point on the image plane. 

Since the state is specified at one of the bound- 
aries and because the system has a forcing input 
G, postulate a solution for costate A as : 


S(x P (r)) = 7, K(x p (r)) = -76? (24) 

Differentiating the expression (23) and substitut- 
ing for G and A from (21) and (22) yields 

A 2 

-a(G - G) = SG - S—[SG + A’] + K (25) 

P 

Separating terms 


A 2 A . a2 

S— 1\ + aG - K = G(S - S 2 —- 


a 


+ 0) (26) 


Since expression (26) must hold for arbitrary 
values of G\ one has that 


.9 = S 2 — - a 


.. ,A 2 

A — S —j- A T oG (28) 

with the boundary condition being given by (24). 

Expression (27) is a Riccati differential equa- 
tion, while (28) is a linear differential equation. 
Once the value of S and K are known as a func- 
tion of .T p , the range can be computed as 

[SG + K]A 

(29) 

The equation (29) is a linear feedback law for 
range £ in terms of the state G and the weighting 
factors. In a following subsection, implementation 
of the range computation algorithm will be given 
in further detail. 

Note that two-point boundary-value problems 
such as the one discussed above can be derived for 
every pair of lines in the registered stereo pair. Ad- 
ditionally, if the weights a, /?, 7 are assigned cer- 
tain statistical significance, the optimization prob- 
lem discussed in the foregoing can be posed as the 
smoothing problem in estimation theory 13 . 


A — SG + A ( 23) 

The first term takes care of the instantaneous de- 
pendence of the costate on the state G', while the 
second term accounts for the external input G. 

In order that expression (23) satisfies the 
boundary condition (20), it is necessary that 


Second-Order Necessary Conditions 

It is known 13 that the sufficient conditions for 
an optimum requires the second variation to be 
strictly positive. This requirement can be met bv 
satisfying the Legendre-Clebsch and Jacobi nec- 
essary conditions in strengthened form. Each of 
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these necessary conditions will be briefly examined 
for the stereo ranging scheme in the following. 

( 1 ) Legendre-Clebsch necessary condition : 
This necessary condition examines the sign 
of the second partial derivative of the variational 
Hamiltonian with respect to the control variable z. 
Thus, a necessary condition for optimality is : 


d 2 H 
dz 2 


> 0 


(30) 


This condition can be found to be satisfied if 


/J>0 (31) 

Since the parameter ft is a user specified 
quantity, the Legendre-Clebsch necessary condi- 
tion does not impose any restrictions on the com- 
putational scheme. 

(2)Jacobi’s necessary condition : 

The stereo ranging process described in the 
foregoing is a linear optimal control problem. 
Hence, Jacobi’s necessary condition will be satis- 
fied with a margin if the differential equations (27) 
and (28) produce finite values and if the variable S 
is positive 13 . Thus, if the range estimates emerging 
from the present computations are finite, Jacobi’s 
necessary condition will be met with a margin. 

In summary, if the parameter (3 is chosen to 
be positive and if the resulting the range estimates 
remain finite, the results would be optimal with 
respect to the specified performance index. 

The Computational Algorithm 

Since the ranging scheme will be tested using dis- 
cretized images, the algorithm implementation will 
be discussed in terms of such images. Consider an 
image consisting of m x n pixel array. Let the hor- 
izontal field-of-view be ^ radians and the vertical 
field-of-view be 6 radians, with / as the camera fo- 
cal length. As the stereo image pair is registered, 
the analysis can be carried out separately for m 
pair of lines. Along any picture line, the distance 
between the centers of any two pixels is given by : 

6 = 2/ta " r/ ' 2 (32) 

n 

Various differential equations involved in the 
stereo algorithm can now be discretized using this 
distance. 

The first step in the computational algorithm 
is the formation of quantities 


and 


, _ 2 (£ 1 -£ 2 ), 

' ' /A* 


( 33 ) 


Gi — (E\ + E 2 )U (34) 


for i — 1,2,3 n, at every pixel along each of the 

m lines of the registered stereo pair. These quan- 
tities are stored for subsequent use. The camera 
focal length / and the stereo baseline Ax are con- 
stants of this process. If the coefficient A(x p ) is 
zero anywhere on the image plane, it can be re- 
placed by a very small value in order to ensure 
complete controllability of the process defined by 

(13). 

Next, the difference equations 


s j+ 1 = Sj + - at 


(35) 


A 2 S 

Aj-j-i = hj + Sj — A — A j + (xGjb (36) 


for j = n, n - 1 , n — 2 , 1 , are propagated from 

the right edge of the image to the left edge using 
the specified conditions on the right edge S n = 
7 , K n = — 7 G n . The results are then stored in 
two arrays. 

As the third step, propagate the differential 
constraint 


Gk+i — Gk — -~—[SkGk + Kk] (37) 

for k = 1 , 2 , ..n, using the condition G\ = G 1 from 
the left boundary of the image to the right edge of 
the image. Store the Gk history in an array. 

The final step is the computation of the range 
using the stored histories S, A\ G using the expres- 
sion 


[SkGk + Kk]Ak 

z k = -jj (33) 

This procedure is repeated for every m hori- 
zontal line in the two images. A flow chart illus- 
trating various steps in the range computations is 
shown in Figure 2 . In the next section, the use of 
this stereo ranging scheme will be illustrated using 
a laboratory image pair. 
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Figure 2: The stereo ranging scheme. 





Figure 4: The laboratory stereo pair: right image. 


Results and Discussions 

The ranging algorithm was tested using a lab- 
oratory stereo pair shown in Figures 3 and 4. 



Figure 3: The laboratory stereo pair: left image. 


The scene consisted of a dark wall in the back- 
ground with a table in the foreground. A soda can 
together with two pencils were then arranged on 
this table. A wire is strung across the two pen- 
cils to explore the feasibility of detecting obstacles 
with small dimensions. A plan view of the various 
objects within the field-of-view is given in Figure 


5. The stereo images were produced using a CCD 
(Charge Coupled Device) camera with a field-of- 
view of about 45 degrees and a focal length of 
about 6 mm. The stereo baseline was 0.1 inches. 
The gray scale consisted of 256 levels correspond 
ing to 8 bit digitization of the CCD image. The 
CCD output was digitized and processed on a SUN 
workstation. 
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Figure 5: A plan view of the various objects in the 
field-of-view. 

Sum and difference between the two images in 
the stereo pair are next formed. These are shown 
in Figures 6 and 7. Range calculations are car- 
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ried out using the data available in these images. 
The weighting factors used in generating the re- 
sults given in this paper are a = 1.0 , (3 = 0.001 , 
7 = lx 10 6 . The solution to the Riccati equation 
(35) and the equation (36) generated using these 
weighting factors are illustrated in figures 8 and 9. 



Figure 6: The laboratory stereo pair: the sum of 
the irradiances. 


Figure 8: Solution of the Riccati equation, S . 




Figure 7: The laboratory stereo pair: the differ- 
ence between the irradiances. 



Figure 9: Solution of the equation (36), I\ . 

These images are scaled so that black corresponds 
to zero and white corresponds to the maximum 
value. Note that the value of S is always greater 
than zero. The results of the range calculation are 
summarized in Table 1. Range estimates are accu- 
rate to within 10% in most cases. To illustrate the 
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Table 1 : Comparison Between Actual Range 

(inches) and the Estimated Range (inches) 


Object 

Actual 

Range 

Estimated 

Range 

No. of 
Points 

Left Pencil 

30.0 

24.13 

911 

Right Pencil 

26.0 

23.37 

1889 

Tape-on-wall 

60.0 

57.82 

639 

Soda can 

28.5 

21.43 

2139 

Bracket 

19.0 

18.55 

982 


nature of range output corresponding to the soda 
can, a histogram showing the number of points at 
one inch intervals is given in Figure 10 . The asym- 
metric nature of the range estimate distribution 
can be seen from this picture. The range com- 



Figure 10: Range histogram for the soda can. 

puted at points of small irradiance differences are 
suceptable to noise. As a result, the ranges corre- 
sponding to such points are dropped from further 
consideration. From Table 1 it may be observed 
that the computed ranges compare favorably with 
the actual ranges. The number of range points 
correponding to various objects show the dense na- 
ture of computed range points. Typically, range is 
computed at about 10 % of the points on the im- 
age plane. To further illustrate the performance 
of the ranging algorithm, the computed ranges are 
displayed in a 3-D format in Figure 11 . Relative 



Figure 11 : Computed ranges from the laboratory 
image pair. 

placement of various range points can be visualized 
from this figure. In Figure 12, the sum of the irra- 
diances produced by the differential constraint (13) 
while computing the ranges is illustrated. Figures 



Figure 12 : Sum of the irradiances predicted by the 
ranging algorithm. 

12 and 6 can be compared to evaluate the effect of 
the weighting parameters a, fl , 7 . Image in Figure 
12 is a smoothed version of image 6 . The weight- 
ing parameters a, /3, 7 should be chosen to retain 
adequate fidelity between Figures 6 and 12 . 

The accuracy of the present range calcula- 
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tion can be improved by employing a finer gray 
scale quantization. Further improvement can be 
achieved by increasing the pixel density in the im- 
age plane. Additionally, one may implement a 
ranging algorithm based on higher-order Fade’ ap- 
proximation of the correspondence equation. Note 
that in the foregoing, the constraint that the com- 
puted range should be greater than zero was not 
explicitly included. Modification of the present for- 
mulation to include this requirement will be a fu- 
ture research item. 

Conclusions 

Vision- based stereo ranging was cast an op- 
timal control problem in this paper. The cor- 
respondence hypothesis was approximated using 
Pade’ approximation and was used as a differen- 
tial constraint in an optimization problem with a 
quadratic performance index. The state variable 
in this formulation was the sum of the image irra- 
diances and the control variable was the range to 
various objects within the field-of-view. 

Since this differential constraint is only an ap- 
proximation and due to the fact that the observed 
irradiances are corrupted by noise, it was satis- 
fied in a minimum error sense by defining a perfor- 
mance index. The performance index consisted of 
a weighted linear combination of the square of the 
irradiance error between the actual irradiance and 
that predicted by the correspondence hypothesis, 
and the range from the image plane. 

The necessary conditions for optimality yield 
linear two-point boundary value problems which 
were solved using a version of backward sweep 
method. Second-order necessary conditions ensur- 
ing the optimality of the range estimates were ex- 
plored. A discrete implementation of the algorithm 
was then discussed. Finally, the ranging algorithm 
performance was demonstrated using a laboratory 
image pair. 

The contributions of the present research are 
the following 

• A new approximation for correpondence hy- 
pothesis in stereo ranging was developed using 
Pade’ expansion. 

• The image-based ranging problem was formu- 
lated as an optimal control problem, obviating 
the need for computing partial derivatives as 
in previous research. 


• It was shown that the optimality of the range 
estimates can be assured under certain mild 
restrictions. 

• The performance of the stereo ranging algo- 
rithm was demonstrated using a laboratory 
stereo pair. 
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6. Guidance Law for Vision- Based Aircraft Maneuvers 


The range data produced by vision-based ranging algorithms are in discrete form. Typically, the 
data is available only at about 10 % of the points. A need exists for deriving guidance laws that 
use the sparse range data to generate a clear trajectory joining the initial and final conditions. An 
exploratory look in this direction is discussed in the following paper. 

In this work, an optimal control problem that maximizes the sum of the distances to various range 
points within the field of view subject to differential constraints arising from the helicopter point-mass 
equations of motion are considered. The performance index additionally includes the square of the 
vehicle acceleration and the flight time, both of which serve to ensure that the resulting trajectories are 
implementable. Feedback linearization of the helicopter dynamics is used to make the present optimal 
control problem amenable to analysis. Necessary conditions for optimality yield a linear two-point 
boundary- value problem which can be solved in feedback form. Inverse transformation of the result 
then produces the desired guidance law. The present analysis resulted in the synthesis of a nonlinear 
guidance law in state feedback form. A paper based on this research was presented at the 1991 AIAA 
Guidance , Navigation, and Control Conference . This paper is given in the following. 
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Abstract 

An optimal guidance scheme for vision-based ob- 
stacle avoidance is developed. The proposed approach 
is useful for automating low-altitude rotorcraft flight. It 
explicitly accounts for the discrete nature of range in- 
formation available from vision- based sensors and uses 
a linear combination of flight time, square of the vehicle 
acceleration and the square of the distance to various 
sensed obstacles as the performance index. A sixth- 
order, three-degree-of-freedom nonlinear point-mass ve- 
hicle model is included in the analysis. Numerical re- 
sults using a sample image sequence is given. 


Introduction 

Passive vision-based obstacle-detection sensors are 
currently being developed by NASA for use with au- 
tomated nap-of-the-earth rotorcraft flight 1 " 6 . The ob- 
jective is to automate the low-altitude rotorcraft flight 
regime using vision-based sensors such as low-light tele- 
vision and infrared imagers. The first aspect of this 
problem, viz., the problem of ranging using motion and 
stereo image sequences is currently being addressed by 
several researchers 2 " 6 . 

The focus of the present paper will be the develop- 
ment of a guidance law for synthesizing flyable trajec- 
tories using the vision-based sensor derived range data. 
Due to the discrete nature of the current imaging sen- 
sors and because of the fact that the range can be com- 
puted only at those regions where a contrast exists, the 
range information is typically available only at about 
10 % of the points in an image. This necessitates the 
development of obstacle-avoidance guidance laws that 
use the discrete range data to synthesize implementable 
trajectories. 

The robotics literature is replete with obstacle- 
avoidance algorithms, most of which have a heuristic 
basis; see Reference 7 for example. Algorithms that 
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use various approximations for obstacles 8 and those 
that require extensive numerical calculations 9 have also 
been discussed in the literature. With the exception of 
the heuristic algorithm discussed in Reference 10, these 
guidance laws are primarily useful for two-dimensional 
maneuvers and assume the availability of the exact ob- 
stacle dimensions and location. 

In the present research, optimal control theory 11 
will be used in conjunction with feedback linearization 
to synthesize a guidance law for obstacle avoidance. 
The research will employ a nonlinear point-mass he- 
licopter model. It will be assumed that the initial and 
final vehicle position and velocity are specified. Using 
a performance index consisting of a linear combination 
of the flight time, square of the rotorcraft acceleration 
magnitude and the sum of the distances to various ob- 
stacles, an optimal control law will be obtained in closed 
form. For this initial study, it will be assumed that the 
obstacles are fixed. Generalization to the case of moving 
obstacles is not difficult, although non-trivial. Further, 
the vehicle as well as the obstacles will be represented 
by points. Using a previously discussed methodology 12 , 
it is possible to tailor the present formulation to include 
the geometric dimensions of the helicopter and the ob- 
stacles. This aspect of the obstacle avoidance guidance 
will be investigated in a future research. Finally, it 
may be noted in passing that the present methodology 
is useful for other vision-based guidance tasks such as 
spacecraft docking and autonomous vehicle guidance. 


Vehicle Model 

The three-degree-of-freedom point-mass model for 
a rotorcraft is given by the following six first-order non- 
linear differential equations : 


, Tsin 0 - D 
V = g sin 7 


X = 


m 

T cos B sin <f> 
mV cos 7 


t _ T cos 0 cos <j> g cos 7 
7 “ mV V~ 

x' = V cos y cos x 


( 1 ) 

( 2 ) 

(3) 

(4) 
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y' = V cos 7 sin x (5) 

z f — V sin 7 (6) 

Here, x is the down range, y the cross range, z alti- 
tude, V airspeed, 7 the flight-path angle, \ heading 
angle, T main-rotor thrust, m the vehicle mass and g 
is the acceleration due to gravity. A prime over a vari- 
able indicates differentiation with respect to time. The 
variable D is the vehicle drag, determined using the at- 
mospheric density p , airspeed V , the reference area s 
and the drag coefficient Cd as: 

D = l -pV 2 sC D (7) 

The drag coefficient for the UH-60 helicopter from Ref- 
erence 13 will be employed in this work. The coordinate 
system used in deriving the vehicle model is illustrated 
in Figure 1. The control variables in these equations 



Figure 1: The coordinate system 

are the rotor thrust orientation angle 0 measured with 
respect a plane normal to the velocity vector, the bank 
angle <f> measured with respect to the local vertical, and 
the main rotor thrust magnitude T. 

If one attempts to derive an obstacle-avoidance 
guidance law using the nonlinear vehicle model de- 
scribed in the foregoing, the solution can only be ob- 
tained using numerical methods. The objective here, 
however, is to obtain a feedback law. To meet this ob- 
jective, the nonlinear system will be first transformed 
to a linear time-invariant form using a coordinate trans- 
formation. This transformation implicitly assumes the 
availability of the vehicle velocity vector components 


from on-board measurements. The guidance problem 
is then solved in the transformed coordinates. Inverse 
transformation of the result then produces the nonlin- 
ear feedback guidance law. 

Differentiating the down-range, cross-range, alti- 
tude equations once with respect to time, the rotorcraft 
equations of motion can be expressed in the following 
form. 


by 


= 3/" = a y , z" = a z (8) 

The right hand sides of these equations are given 


a s = V* cos 7 cos x — X*V cos 7 sin \ 

— y'V sin 7 cos x (9) 

a y = V* cos 7 sin \ + x*V cos 7 cos \ 

-7 / Ksin7sinx (10) 


a z = V* sin 7 + 7' V cos 7 (11) 

The variables V f ,y\x' ma y be eliminated from equa- 
tions (9) - (11) using expressions (1) - (3). 

The vehicle acceleration components a x , a y , a z will 
be treated as the new control variables in the guidance 
problem and will be termed the pseudocontrol vari- 
ables in the ensuing. Note that the system dynamics is 
linear in terms of the pseudo-control variables. If these 
variables were known, together with the vehicle velocity 
vector components x', j/, z ’ , the actual control variables 
can be computed as: 


^tan- 1 K cos *7 Q * sin * 

l r 


9 — tan" 1 CQS< ^(^ + ^2) -f A3 
r 


T = 


mT 


where, 


cos 9 cos <j> 

V = y/x' 2 + y' 2 + ~ 2 
7 = sin -1 (z'/V), x = tan ~ l (y'/x') 


( 12 ) 

(13) 

(14) 

(15) 

(16) 


Ni = a x cos 7 cos x + o, y cos 7 sin x (17) 


N2 = (a z + g)siny } N 3 — Dcos<j>/m (18) 


r = (g + a z ) cos 7 — a x sin 7 cos X — a y sin 7 sin x (19) 

The vehicle model is now in the form of a linear 
time-invariant system. If the obstacle avoidance prob- 
lem is formulated in transformed coordinates, it can be 
solved in closed form. 
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Guidance Law Development 

Given the current vehicle position and the desired 
final vehicle position, the objective is to produce a tra- 
jectory that maximizes the sum of the square of the dis- 
tances from the vehicle to the various obstacles located 
at the points x,* , t/» , Z { , * = l,2,...n. Further, to enhance 
the trajectory masking, the vehicle is required to stay 
as close as possible to a specified altitude z r while using 
the least square of vehicle acceleration magnitude. It 
is sometimes desirable to execute the mission in mini- 
mum flight time. Thus, a composite performance index 
of the following form may be specified. 


min / / {c + 
> a v> a * Jo { 


( 1-0 


a(a x 7 + a v 7 + a 2 7 ) 


+K Z ~ z rf ~ 0 5 ^[(* “ 


i = 1 

+(y - in) 2 + (*- z >) 2 ] 


dt 


( 20 ) 


= (25) 

i= 1 


= (1 - C) 


0^2( z - *i) - K Z - Zr) 


<■ i= 1 


(26) 


A X ' Ay' \ z ' 

ax = ^7) ’ = = ^~-7) 

The strict upper bound on the parameter £ in equation 
(21) is required to ensure that the control variables re- 
main finite. 

Since all the states are specified at the initial and 
final time, no boundary conditions can be specified on 
the costates. The state-costate system for this problem 
is linear. Consequently, the closed-form solution can 
be obtained by proceeding as follows. Differentiate the 
A'r', A'y/, y equations (23) with respect to time and 
substitute for A' x , A' y , A'* from the expressions (24) - 
(26). This yields : 


The initial vehicle position components xq, t/o, 2o, the 
velocity components x'o, t/ 0 > z *o , the obstacle position 
components x,, y,, z iy i = l,2,....n, the final vehicle 
position components xy, yy, Zy , and the velocity com- 
ponents x'y, j/y, z* y are assumed to be specified. The 
positive parameters f , a, fi, ft in the performance index 
establish the relative weighting of flight time, vehicle 
acceleration, altitude deviation from the reference alti- 
tude and the distance to various obstacles. For reasons 
that will be made clear in the ensuing, the parameters 
C, ft and fi are required to satisfy the constraints 

0<C<1, nft>t* ( 21 ) 

The variational Hamiltonian 11 for this problem is 
defined as 


H = C + 


(i-C) 

2 



+ o u 


+ a, 7 ) 


+n(z - z T ) 7 - 0 - x,) 2 

i = 1 

+(y - y;) 2 + {z- 2 .) 2 )| + K'Qx 

4-Ay'ay + A z >a z *+■ A r x / H~ A y y + A Z z (22) 
The necessary conditions for optimality are 

A' r / = -A,, AV = -A y , y z > = -\ z (23) 


A% = -(!-() 


H 

ft- 

2; 

vy 

1 

(28) 

1 

% 

<c; 

1 

(29) 

i = l 


n 

~ *•') “ V{Z - Zr) 

1 = 1 

(30) 


Next, substitute for a x ,a y ,a z in terms of the 
costates A r s A y /,A 2 ' in the x" , y", z" equations to 
yield: 


—A X ' 

«( i-O’ 


— Ay' 


a(l-C) 


(31) 


Differentiating the equations (31) twice with re- 
spect to time and substituting for A" r ', A ;/ y /, A " z > from 
(28) - (30) results in three fourth-order linear differen- 
tial equations of the form 


//// 0 r 1 

i = l 

y"" = ^i n y - J2 w] 


]tUL z _ty z . + fL 

a ' a 

» = i 


a 


z r 


(32) 

(33) 

(34) 


v* = (i - 

i — 1 


(24) 


Denoting a x - 
tions (32) and (33) are 


, the eigenvalues of the equa- 
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\ J n ^ a ^ > the eigenvalues for the equation (34) are given 
by : 


*i, -<rijau -M (35) 

Note that j = y/^1. Similarly, denoting <T 2 = 


^ 2 , “< 72 , . 7 ^ 2 , - J <^2 


(36) 

The constraint on the product n/3 as in equation 
(21) is essential for the obtaining the real eigenvalues 
in (36). 

Now, from the theory of linear differential 
equations 14 , the general solutions to the necessary con- 
ditions (32) - (34) can be obtained as: 


x(f) = A x e* xi + B x e ° lt + C x sin a\t 

+D X cos ait + s x (37) 

y(t) = A y e 0lt + B y e~° xt + C y sin ait 

-f-D y cos ait + s y (38) 

z(t) = Ate** 1 -f Bte"** 1 + C z sin a 2 t 

+D Z cos a 2 t + s z (39) 


where 


«=1 1=1 


1=1 

1 


nfi- H 


• 1 = 1 


(40) 

(41) 


Ax , Ay , Az , B X , By , Bz , C X , Cy , Cz , D x , Dy , Dz 
are the constants determined using the specified bound- 
ary conditions. It may be verified that the costates in 
this problem are given by: 


A*' = a« - l)crj : 


A x e ait + Bxe- 0 '* 


—C x sin (7 — D x cos a\t 


(42) 


= a(C - l)<rC 


A y e °' t + B y e~° lt 


—Cy sin (T \t — D y cos (?\t 


(43) 


A,* = o(C - l)^ 2 


A 2 e° 3t + B 2 e _<,3 ‘ 


— C 2 sin a 2 t — D z cos < 72 t 

\ x = a(l - CVi 3 ]Axe a - B x e~ a>t 
— C x cos ait + D x sin ait 

\y=a{\-Qa l 3 

—Cy cos (r \t + Dy sin <r\t 
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(45) 


Aye ait -Bye-' 1 * 


(46) 


A 2 = a(l - C)o - 2 3 


Aze Clt - B 2 e-° 


—Cz cos a 2 t + D 2 sin ait 


(47) 


Next, applying the given boundary conditions 
along the x,y, z axes leads to the matrix equations 


(48) 
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(52) 

(53) 


52 = sin<7 2 f/, C 2 = cos a 2 t / (54) 

The matrix equations (48) - (50) can be solved to 
yield the 12 arbitrary constants in the problem if the 
final time tj were known. 

The final time can be determined by invoking a 
constant of motion in this problem. Since the varia- 
tional Hamiltonian is autonomous, and since the time 
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weight parameter £ is constrained to be greater than 
zero, the final time may be found using the transversal- 
ity condition: 


H(t) = H(t f ) = 0 (55) 

Satisfaction of this condition would involve a one- 
dimensional search wherein one would assume a final 
time and determine the value of the variational Hamil- 
tonian by first finding the arbitrary constants and then 
substituting for the states, costates and the control vari- 
ables. If this results in the satisfaction of equation (55), 
then the optimal final time has been found. Otherwise, 
update the assumed final time and repeat these calcu- 
lations. One-dimensional search techniques such as the 
method of bisections can be used to rapidly determine 
this quantity. 

Note that at this stage, the solution is in open- 
loop form. Since the position of the obstacles and the 
vehicle states may not be known exactly, a closed-loop 
solution is desirable for on-board implementation. To 
this end, assume that the current time is the initial 
time, f = 0. The current vehicle states then become 
the initial conditions and the final time t j becomes the 
time-to-go, denoted by the variable t go . The current 
commanded vehicle accelerations are then given by the 
expressions 


a x - a i 2 
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k 8 = 


<7 2 (sinh aitg 0 - sin aiig 0 ) 
cos aitgo cosh aitgo 1 


(69) 


The only unknown parameter in these equations is the 
time-to-go. This quantity can be computed using a 
one- dimensional numerical search just as in the com- 
putation of */ discussed previously. 

Once the variables a Xl a y , a z are calculated, the real 
control variables of the aircraft can be recovered using 
the transformations given in equations (12) -(14). Fur- 
ther details of the guidance law implementation will be 
described in the next section. 

Note that the controls emerging from the present 
development satisfy the strengthened Legendre-Clebsch 
necessary condition. Verification of the Jacobi’s neces- 
sary condition will be a future research item. 


The unknowns in the expression (56) - (58) are the 
arbitrary constants D Xy D y) D z . These may be found 
using the matrix equations (48) - (50). The resulting 
closed-loop guidance law will be of the form 
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where 
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(62) 


Algorithm Implementation and 
Evaluation 

The obstacle-avoidance guidance law was imple- 
mented on a nonlinear point-mass simulation of the ro- 
torcraft. The salient steps in implementing the guid- 
ance law are: at each guidance interval, 

1. Select the acceleration weight a, the obstacle- 
avoidance weight /?, the altitude-deviation con- 
straint weight fi and the time weight C 

2. Compute the eigenvalue magnitude <7i, <7 2 , and 
Sx^Sy ] s 2 using the given formulae. 

3. Using the current states and the final states, com- 
pute the time-to-go using the transversality condi- 
tion H(t 0 ) = 0. 


ki 


<7! 2 (cos aitgo — cosh aitgo) 
cos aitgo cosh a x t go — 1 


(63) 
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4. Compute the feedback gains k\...kg and calculate 
the desired vehicle accelerations a x ,a yi a z in feed- 
back form. 



5. Using the commanded vehicle acceleration, com- 
pute the required rotorcraft main rotor thrust T, 
thrust orientation angle 6 and the bank angle <f>. 

In order to evaluate the performance of the guid- 
ance law, trajectories were generated using a range map 
obtained using the algorithm discussed in Reference 2. 
Figure 2 shows the sample scene obtained from a cam- 
era fixed on the nose of a helicopter flying at an altitude 


Time Sync. Data 

Helicopter Nose Boom 



Figure 2: A Sample Image 


of about 7 feet. Besides other things, the scene shows a 
runway, and four vehicles parked on the two sides. The 
instrumented nose boom of the helicopter is clearly visi- 
ble at the top right hand part of the image. Parameters 
identifying this frame are displayed in a sub-picture on 
the top left hand side of the image. 

This image consisted of a 512 x 512 pixel array, 
with 8 bit gray-scale digitization. The range data was 
obtained by processing several of these airborne cam- 
era images using the vision-algorithm discussed in Ref- 
erence 2. This algorithm produced range data at the 
points denoted by small white squares in Figure 3. The 
range data was available at 293 points on the image 
plane. This represents about 0.1 % of the possible 
262,144 points. The range data is recovered only at 
those points where a significant image intensity change 
is detected between frames. However, note that the 
algorithm determines the range to the four vehicles 
parked on the sides of the runway. 

For the purposes of illustration, the mission to be 
flown was assumed to be that of a helicopter transi- 
tioning from forward flight to a landing mode with zero 
vertical velocity and small horizontal velocity. The tar- 
get touch-down is chosen to be at a point between the 



Figure 3: Available Range Data 


stationary vehicles 1 and 2. The boundary conditions 
corresponding to this trajectory were as follows. The 
initial velocity vector was [0.0, —38.0, 0.0] T feet/s, and 
the vehicle position vector was [532.0, 590.0, 7.0] r feet. 
The specified final velocity and position vectors were 
[1.0, 0.0, 0.0] r feet/s , and of [740.0, 0.0, 0.0] T feet 
respectively. The weighting factors in the performance 
index were chosen as ( = 0.96, a = 0.5, /i = 0.0,/? = 
10“ 06 . The flight time for this trajectory turns out to 
be 20 seconds. 

The evolution of the rotorcraft altitude as a func- 
tion of down range is shown in Figure 4. Although the 
present analysis considers only point obstacles, these 
obstacles have been denoted by small circles for the sake 
of clarity. The horizontal projection of the trajectory 
is shown in Figure 5. In order to illustrate the obstacle 
avoidance charecteristics of the trajectory, the distance 
to the nearest obstacle at each time instant along the 
vehicle trajectory is plotted in Figure 6. In the present 
example, the minimum distance was about 8.25 feet. 
Note that it may be possible to further increase this 
clearance by increasing the obstacle avoidance weight 
parameter /?. The rotorcraft airspeed along the trajec- 
tory is shown in Figure 7. The vehicle initially accel- 
erates through the turn before decelerating to satisfy 
the specified terminal velocity. The temporal evolution 
of the rotorcraft heading angle is illustrated in Figure 
8. The guidance law commands a large heading angle 
correction towards the end to meet the specified final 
boundary condition on the airspeed. Throughout the 
maneuver, the vehicle load factor was within 1.02, and 
the maximum magnitude of the rotor thrust orientation 
angle 0 was less than about 12 degrees. The bank angle 
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4> was well within 15 degrees as may be observed from 
Figure 9. To get a feel for the nature of trajectories 
produced by the guidance law, the rotorcraft trajec- 
tory is next superimposed on the sample image using a 
perspective projection in Figure 10. The point marked 
with an inverted *T* denotes the current position of the 
helicopter. The vertical lines are used to indicate the 
helicopter altitude above the runway surface along the 
trajectory. The obstacle-avoidance characteristics are 
apparent from this figure. 


Conclusions 

Development of an obstacle-avoidance guidance 
law that uses vision-based range data was presented. 
This analysis used a sixth-order nonlinear point-mass 
model of the helicopter together with a linear combina- 
tion of flight time, square of the rotorcraft acceleration 
magnitude and the square of the distance to various ob- 
stacles as the performance index. The obstacles were 
represented as points. 

The rotorcraft model was first transformed into lin- 
ear, time-invariant form using a coordinate transforma- 
tion. The guidance problem was then solved using the 
transformed model. The resulting guidance law is in lin- 
ear feedback form. Inverse transformation of the guid- 
ance law yields the vehicle guidance commands. The 
performance of the guidance law was illustrated using 
a realistic vision-derived range data. 

The present guidance law is useful for other vision- 
based guidance tasks such as spacecraft docking and 
autonomous vehicle guidance. 


Figure 9: Optimal Trajectory: Rotorcraft Bank Angle 
Vs Time 



Figure 10: Perspective Projection of the Optimal Tra- 
jectory on the Image Plane 
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7. Conclusions 


This report outlined the research carried out on the development of field-based ranging algorithms 
for rotorcraft nap-of-the-earth flight. Based on various publications that have resulted from this work, 
the contributions may be summarized as follows : 

• A field-based method for ranging using motion image sequences was developed by combining the 
Horn-Schunk image constraint equation with expressions for incremental perspective projection 
and an irradiance tracking filter. This scheme was tested using a simulated image sequence 
and included the translational camera motion. The spatial-temporal sampling requirements to 
obtain a specified ranging accuracy were examined. 

• The field-based ranging algorithm was next generalized to include both motion and stereo image 
sequences by replacing the Horn-Schunk image constraint equation with a multi-dimensional 
Taylor series approximation for the correspondence hypothesis. This step produces a set of 
ranging equations, together with expressions that predict the error involved in the Taylor series 
approximation. Several orders of the ranging algorithm was tested using laboratory image se- 
quences. The lowest-order approximation was found to be adequate in most image sequences 
collected in the laboratory. 

• The ranging algorithms require estimates of the spatial partial derivatives of the image irradi- 
ances. A method for estimating partial derivatives by product factorization of the images was 
developed. This method converts the partial derivative estimation problem into a set of linear 
lumped-parameter estimation problems, This method was tested on several laboratory images 
and found to produce excellent partial derivative estimates. The factorization approach for par- 
tial derivative estimation was next used in conjunction with the ranging equation to yield a fast 
stereo ranging algorithm. 

• By defining various coordinate systems and the incremental transformations, the image-based 
ranging algorithm was extended to include rotational and translational motion of the rotorcraft 
and the cameras. Both ranging equation as well as the range error equation were developed. 
It was shown that ranging can be accomplished if the incremental translation and incremental 
rotation angle of the camera are known. 

• The need for the computation of the partial derivative before ranging prompted research on 
methods that do not require partial derivative estimates. Specifically, the stereo ranging problem 
was examined. First, it was shown that a Pade approximation can be used to approximate 
the correspondence hypothesis. Depending on the nature of the image sequence, it might be 
useful to employ different orders of Pade approximation. For the stereo problem, the first-order 
approximation produces a varying coefficient first-order linear ordinary differential equation. 
This ordinary differential equation turns out to be identical to that synthesized by combining 
the backward and forward Taylor series approximation of the correspondence equation. Instead 
of satisfying this differential equation by first computing derivatives, an optimization problem 
was posed whereby the irradiance predicated by the Pade approximation is compared with that 
from the actual irradiance. Ranges that minimize the integral of the irradiance error along each 
image line is then found using the necessary conditions for optimality. 

• Finally, research on using the discrete vision-based range data for optimal vehicle guidance was 
initiated. The problem of optimally navigating a rotorcraft through a field of point obstacles was 
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considered. This research used a point-mass vehicle model and a quadratic criterion. Necessary 
conditions for optimality are then used to obtain a nonlinear feedback guidance law. 

A few papers were presented at various national conferences based on the present research. These 
are listed below in the order in which they appeared. 

[1] Menon, P. K. A., and Sridhar, B., “Passive Navigation Using Image Irradiance Tracking,” AIAA 
Guidance, Navigation, and Control Conference , August 14-16, 1989, Boston, MA. Parts of this work 
was also presented at NASA Vision Science and Technology Workshop, Nov. 30 - Dec. 2, 1988, NASA 
Ames Research Center, Moffett Field, CA; and AIAA Houston Chapter Invitational Conference on 
Guidance and Control, NASA Johnson Space Research Center, February 12, 1990. 


[2] Menon, P. K. A., and Sridhar, B., “Image-Based Range Determination,” AIAA Guidance, Naviga- 
tion, and Control Conference , August 14-16, 1990, Portland, OR; Also being revised for the Journal 
of Guidance, Control , and Dynamics. Parts of this work was also presented at NASA Workshop on 
Vision-Based Rotorcraft Navigation, September 19, 1990. 

[3] Menon, P. R. A., Chatterji, G. B., and Sridhar, B,, “A Fast Algorithm for Image-Based Ranging,” 
SPIE International Symposium on Optical Engineering and Photonics in Aerospace Sensing , April 
1-5, 1991, Orlando, FL. 


[4] Menon, P. K. A., Chatterji, G. B., and Sridhar, B., “Passive Obstacle Location for Rotorcraft 
Guidance,” AIAA Guidance , Navigation , and Control Conference , August 12-14, 1991, New Orleans 
LA. 

[5] Menon, P. K. A., Chatterji, G. B., and Sridhar, B., “Electro- Optical Navigation for Aircraft,” 
Submitted for consideration in the IEEE Transactions on Aerospace and Electronic Systems. 

[6] Menon, P. K. A., Chatterji, G. B., and Sridhar, B., “Vision-Based Optimal Obstacle- Avoidance 
Guidance for Rotorcraft,” AIAA Guidance , Navigation, and Control Conference , August 12-14, 1991, 
New Orleans, LA. 

[7] Menon, P. K. A., Sridhar, B,, and Chatterji, G. B., “Vision-Based Ranging as an Optimal Control 
Problem,” Paper communicated to AIAA Guidance , Navigation, and Control Conference , August 10 
- 12, 1992, Hilton Head, SC. 


With this background, the following items are suggested as promising future research directions. 

1. The ranging algorithms discussed in this report use a pair of images. It is sometimes desirable 
to carry out ranging using several images simultaneously. Indeed, such an approach may yield 
more accurate range estimates. The field-based ranging algorithms need to be reformulated to 
simultaneously handle mutiple images. 

2. The present ranging algorithms uses just the scene irradiances and perspective projection ge- 
ometry to construct range to various objects within the field-of-view. This process completely 
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ignores the relationship between surface brightness, direction of illumination and the surface 
orientation. Invoking appropriate reflectance models such as the Lambert/an surface model in 
conjunction with the ranging equations developed under the present research may yield more 
consistent range estimates. Such an approach would directly produce 3-D surface descriptions. 
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